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1.  Project  Summary 

The  flame  type  of  most  practical  combustion  devices  is  the  diffusion  flame.  These 
flames  are  important  in  the  interaction  of  heat  and  mass  transfer  with  chemical  reactions 
in  ram  jets,  jet  turbines  and  commercial  burners.  Three-dimensional  models  that  couple 
the  effects  of  fluid  flow  with  detailed  chemical  reactions  are  as  yet  computationally  infea¬ 
sible.  We  can,  however,  obtain  important  information  in  practical  systems  by  considering 
two-dimensional  configurations.  In  particular,  in  this  report  we  focus  our  attention  on 
axisymmetric  laminar  and  turbulent  diffusion  flames  in  which  a  cylindrical  fuel  stream  is 
surrounded  by  a  coflowing  oxidizer  jet.  In  this  configuration  we  can  study  the  interaction 
of  fluid  flow  with  chemical  reactions  while  obtaining  a  computationally  feasible  problem. 
The  work  centers  on  the  development  and  application  of  accurate  and  efficient  computa¬ 
tional  methods  for  the  solution  of  the  two-dimensional  nonlinear  boundary  value  problems 
describing  the  reacting  systems.  In  particular,  our  goals  involve  the  generalization  of 
our  one-dimensional  fluid  chemistry  model  to  two  dimensions.  We  also  focus  on  the  use 
of  two-dimensional  flame  sheet  models  as  starting  estimates  for  the  nonlinear  equation 
solver.  The  highly  nonlinear  finite  rate  chemistry-fluid  equations  are  solved  by  a  combina¬ 
tion  of  time  integration  and  a  modified  Newton  method.  The  final  object  of  the  research 
is  the  incorporation  of  a  k  -  c  turbulence  model  into  the  finite  rate  chemistry  formulation. 
Both  confined  and  unconfined  methane-air  flames  have  been  studied.  The  results  of  the 
research  are  applicable  to  problems  in  1)  turbulent  reacting  flows,  2)  engine  efficiency,  3) 
commercial  power  generation  units,  and  4)  pollutant  formation. 

2.  Identification  and  Significance  of  the  Problem 

The  ability  to  predict  the  coupled  effects  of  complex  transport  phenomena  with  de¬ 
tailed  chemical  kinetics  is  critical  in  an  understanding  of  turbulent  reacting  flows,  in  im¬ 
proving  engine  efficiency  and  in  the  study  of  pollutant  formation.  Since  three-dimensional 
models  combining  both  fluid  dynamical  effects  with  finite  rate  chemistry  are  as  yet  com¬ 
putationally  infeasible,  the  modeling  of  chemically  reacting  flows  generally  proceeds  along 
two  independent  paths.  In  one  case  chemistry  is  given  priority  over  fluid  mechanical  ef¬ 
fects  and  these  models  are  used  to  assess  the  important  elementary  reaction  paths,  for 
example,  in  hydrocarbon  fuels.  In  the  other  case,  multidimensional  fluid  dynamical  effects 
are  emphasized  with  chemistry  receiving  little  or  no  priority.  The  goal  in  reacting  flow 
computations,  however,  is  to  be  able  to  combine  the  effects  of  detailed  chemistry  with 
complex  fluid  mechanics. 

3.  Phase  II  Technical  Objectives 

We  have  investigated  the  coupled  effects  of  fluid  flow  with  finite  rate  kinetic  models  in 
two-dimensional  axisymmetric  diffusion  flames.  The  work  centered  on  the  developement 
and  application  of  accurate  and  efficient  computational  methods  for  the  solution  of  the 
two-dimensional  nonlinear  boundary  value  problems  describing  the  reacting  systems.  In 
particular,  our  goals  involved  the  generalization  of  our  one-dimensional  fluid-chemistry 
model  to  two  dimensions.  We  considered  a  stream  function-vorticity  formulation  with 
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detailed  kinetic  theory  transport  and  finite  rate  chemistry.  The  Phase  II  research  plan 
focused  on  the  use  of  two-dimensional  flame  sheet  models  as  starting  estimates  for  the 
nonlinear  equation  solver.  The  highly  nonlinear  finite  rate  chemistry-fluid  equations  are 
solved  by  a  combination  of  time  integration  and  a  modified  Newton  method.  The  final 
objective  of  the  research  was  the  incorporation  of  a  k  -  e  turbulence  model  into  the  finite 
rate  chemistry  formulation.  The  two-dimensional  model  and  solution  procedure  represent 
a  generalization  of  our  Phase  I  initiative.  The  work  proposed  in  Phase  II  has  allowed  the 
investigation  of  the  effects  of  detailed  chemistry  and  transport  on  the  structure  and  extinc¬ 
tion  of  hydrocarbon  flames.  The  researcher  now  has  the  ability  to  predict  the  temperature, 
major  and  minor  species  as  well  as  the  velocities  in  axisymmetric  diffusion  flames.  This 
level  of  predictive  capability  is  useful  to  both  the  Federal  Government  and  the  commer¬ 
cial  sector  in  l)  determining  flame  extinction  limits,  in  2)  predicting  the  effect  of  higher 
hydrocarbons  in  the  onset  of  soot  formation  and  in  3)  an  understanding  of  flame  shape, 
height,  and  width. 

4.  Phase  II  Work  Plan 

The  modeling  of  axisymmetric  diffusion  flames  can  be  reduced  to  the  solution  of  a  set 
of  coupled  nonlinear  boundary  value  problems.  In  these  problems  we  desire  solution  profiles 
to  as  many  as  several  dozen  species  concentrations  in  addition  to  the  temperature  and  the 
velocity  fields.  Although  axisymmetric  flames  are  important  in  combustion  applications, 
they  have  received  relatively  little  attention  in  theoretical  flame  studies.  Part  of  this 
neglect  is  due  to  the  two-dimensional  nature  of  the  problem  coupled  with  the  complexities 
associated  with  the  combined  effects  of  transport  phenomena  and  chemical  processes.  In 
the  axisymmetric  diffusion  flame  we  consider,  a  fuel  jet  discharges  into  a  laminar  air  stream. 
The  flames  can  be  either  confined  or  unconfined.  In  both  cases  the  tubes  through  which 
the  fuel  and  the  oxidizer  flow  are  concentric  and  have  radii  Rj  and  Rq,  respectively.  The 
two  gases  make  contact  at  the  outlet  of  the  inner  tube  and  a  flame  that  resembles  a  candle 
is  produced. 

Our  work  plan  is  divided  into  five  sections.  In  the  first  part  we  consider  the  formu¬ 
lation  of  the  governing  axisymmetric  equations.  We  consider  a  stream  function-vorticity 
formulation.  We  next  incorporate  detailed  transport  phenomena  with  finite  rate  chemistry 
into  the  model.  This  involves  the  use  of  kinetic  theory  expressions  for  the  diffusion  coeffi¬ 
cients,  the  viscosity  and  the  thermal  conductivity.  The  reaction  mechanism  is  be  inputed 
via  CHEMKIN,  the  chemical  kinetics  package  developed  at  Sandia  National  Laboratories 
in  Livermore,  CA  (l].  In  the  third  stage,  multi-dimensional  flame  sheet  models  are  used 
to  provide  starting  estimates  for  the  nonlinear  equation  solver.  In  part  four  we  employ  a 
combination  of  time-stepping  and  a  modified  Newton  method  in  solving  the  large  system  of 
highly  nonlinear  equations.  The  Newton  equations  are  be  solved  by  a  generalized  minimum 
residual  method  for  nonsymmetric  systems  of  linear  equations.  Finally,  we  incorporate  a 
k  —  e  turbulence  model  into  our  fluid  dynamic-thermochemistry  model. 

4.1.  Laminar  Flames 

4.1.1.  Problem  Formulation 

Conclusions  derived  from  studies  of  laminar  flames  are  important  in  characterizing  the 
combustion  processes  occurring  in  turbulent  flames,  in  improving  engine  efficiency  and  in 
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understanding  the  formation  of  combustion  based  pollutants.  By  studying  laminar  flames, 
we  can  identify  the  important  reactions  controlling  extinction  and  we  can  identify  the 
important  species  involved  in  pollutant  formation  while  providing  information  on  the  fluid 
mechanics  of  the  flames.  One  of  the  simplest  two-dimensional  flame  configurations  with 
practical  importance  is  the  axisymmetric  diffusion  flame.  Although  axisymmetric  flames 
are  important  in  combustion  applications,  they  have  received  relatively  little  attention  in 
theoretical  flame  studies.  Part  of  this  neglect  is  due  to  the  two-dimensional  nature  of  the 
problem  coupled  with  the  complexities  associated  with  the  combined  effects  of  transport 
phenomena  and  chemical  processes.  In  this  report  we  consider  two  axisymmetric  diffusion 
flame  configurations-an  unconfined  and  a  confined  flame.  In  the  unconfined  case  a  fuel  jet 
discharges  into  a  laminar  air  stream.  The  tubes  through  which  the  fuel  and  the  oxidizer 
flow  are  concentric  and  have  radii  Rj  and  Rq ,  respectively.  The  two  gases  make  contact 
at  the  outlet  of  the  inner  tube  and  a  flame  that  resembles  a  candle  results  (see  Figure 
1).  In  the  confined  problem  a  cylindrical  shield  surrounds  the  fuel  and  oxidizer  tubes  (see 
Figure  2).  Our  model  of  axisymmetric  laminar  diffusion  flames  considers  the  full  set  of 
two-dimensional  governing  equations.  In  primitive  variables  (r  and  z  denote  the  radial  and 
axial  coordinates,  respectively)  the  governing  equations  can  be  written  in  the  form 
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The  system  is  closed  with  appropriate  boundary  conditions  on  each  side  of  the  com¬ 
putational  domain.  For  the  unconfined  flame  we  have 


Axis  of  Symmetry  (r  =  0): 


Exit  ( z  -*  00): 


dp  dva  dYk  dT  n 

dr  ~Vr  ~  dr  ~  dr  ~  dr  ~  ' 


k  =  1, 2, ...  ,/if, 


(4.7) 


dvr  dvg 
dz  dz 


?Yk 

dz 


dT 

dz 


0,  *  =  1,2 . K, 


(4.8) 


Inlet  (a  =  0): 
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Outer  Zone  (r  =  Rq): 
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Identical  boundary  conditions  are  employed  in  the  confined  flame  except  for  the  outer 
zone.  On  the  shield’s  surface  we  have 


Outer  Zone  (r  =  Rq)- 
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The  subscripts  I  and  O  refer  to  the  inner  jet  and  the  outer  jet,  respectively,  and  pi,po>  vhvO> 
YkfiY^o^Tf, Tq  and  Twau  are  specified  quantities. 

We  can  reduce  the  size  of  the  system  to  be  solved  by  introducing  the  vorticity  and 
the  stream  function  (2).  The  vorticity  is  a  measure  of  the  counterclockwise  rotation  in  the 
flow.  In  particular,  formulation  of  the  vorticity  transport  equation  serves  to  eliminate  the 
pressure  as  one  of  the  dependent  variables.  We  define  the  vorticity  such  that 


dvr  dvg 
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The  stream  function  ^  is  used  to  replace  the  radial  and  axial  components  of  the  velocity 
vector  by  a  single  function.  It  is  defined  in  such  a  way  that  the  continuity  equation  is 
satisfied  identically.  We  have 
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With  the  definitions  in  (4.13)-(4.15),  the  governing  equations  become 
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where  the  components  of  the  iso  operator  are  given  by  ( dpdz ,  —dpdr).  The  boundary 
conditions  for  the  unconfined  flame  in  the  stream  function-vorticity  formulation  are  written 
in  the  form 


Axis  of  Symmetry  (r  =  0): 


Exit  ( z  —*  oo): 


.  dYk 


dT 

dr 


0,  k  =  1,2, 


K, 


(4.20) 


_dui  _  dYk  dT  _ 
dz  dz  dz  dz 


K, 


(4.21) 


Inlet  ( z  =  0): 
r  <  kj 


if)  =  p/u/r2/2, 


R[  <  r  <  Rq 


u)  —  0, 

n  =  n,,  fc  =  i,2,. . . ,k, 

T  =  TIt 

i>  =  \piviR2i  +  \poVoir2  -  R2i)> 

u)  —  0, 

Yk  =  Yko,  k  =  1,2,.  ..,K, 

T  =  T0 , 


(4.22) 


(4.23) 
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Outer  Zone  (r  =  Rq)'- 


=  \pivir /  +  \povo{ro  -  r2i)> 


du 

!h 


=  0, 


dYk 

dr 


=  0,  k  = 


1,2,...,  K, 


(4.24) 


dT 

dr 


=  0. 


The  boundary  conditions  on  the  shield  must  again  be  modified  for  the  confined  flame.  We 
have 


Outer  Zone  (r  =  i?o): 

'P  =  \piviR]  +  ^PovoiRo  ~  r2i)> 

3  07  -  0/-1  _  o'/- 1 

*OP/  (*0  -  R 7_!)2  2  ’ 

"dT=0’  *  =  l,2,...,if,  (4.25) 

r  =  Twau, 

where  quantities  with  the  subscript  I  —  1  in  the  vorticity  boundary  condition  are  evaluated 
at  the  grid  points  next  to  the  shield. 

In  addition  to  the  variables  already  defined,  T,  denotes  the  temperature;  Y^,  the  mass 
fraction  of  the  k ^  species;  p,  the  pressure;  vr  and  vz ,  the  velocities  of  the  fluid  mixture 
in  the  radial  and  axial  directions,  respectively;  p,  the  mass  density;  IV*,  the  molecular 
weight  of  the  k ^  species;  W,  the  mean  molecular  weight  of  the  mixture;  R,  the  universal 
gas  constant;  A,  the  thermal  conductivity  of  the  mixture;  ep,  the  constant  pressure  heat 

capacity  of  the  mixture;  Cp*,  the  constant  pressure  heat  capacity  of  the  k ^  species;  ti>*,  the 
molar  rate  of  production  of  the  species  per  unit  volume;  h*,  the  specific  enthalpy  of  the 
k ^  species;  g ,  the  gravitational  constant;  n  the  viscosity  of  the  mixture  and  and  V*  , 

the  diffusion  velocities  of  the  k^  species  in  the  radial  and  axial  directions,  respectively. 
Specifically,  we  write 


i 


Vkr=-(l/Xk)Dk^,  *=1,2 . K,  (4.26) 

vk,  =  -{UXk)Dt^,  *=1,2 . K.  (4.27) 
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where  Xjt  is  the  mole  fraction  of  the  k ^  species  and  D h  is  related  to  the  binary  diffusion 
coefficients  through  the  relation  (see,  e.g .,  (3|) 


Dk=  ~  V>1  . 


(4.28) 


4.1.2.  Transport  and  Chemistry  Model 

The  transport  of  momentum,  species  mass  fractions  and  energy  in  a  chemically  react¬ 
ing  flow  requires  the  evaluation  of  the  mixture  viscosity,  the  thermal  conductivity  of  the 
mixture  and  the  diffusion  coefficients  of  the  species.  Accurate  and  efficient  evaluation  of 
these  quantities  is  an  important  aspect  of  combustion  calculations.  We  have  modeled  the 
given  transport  coefficients  by  the  following  procedures. 

Diffusion  Velocity 

The  binary  diffusion  coefficients  Djk  can  be  written  as  functions  of  the  temperature 
and  the  pressure  [4].  We  have 


D’k  ~  16 


3  y/^k%ryM~k 


(4.29) 


In  the  above  formula  My denotes  the  reduced  mass  (molecular)  for  the  (j,k)  species  pair 

2  MjMft 


(4.30) 


Oji,.  is  the  averaged  collision  diameter,  is  a  collision  integral  and  kg  denotes  Boltz¬ 

mann's  constant.  A  similar  expression  can  be  written  for  Vj^. 

Viscosity 

The  single  component  viscosities  can  be  written  in  the  form  (4). 


5  \J if Mf/fkgT 

"*  =  16  „„2n(2,J).  > 


(4.31) 


where  M*  is  the  mass  of  the  molecule,  is  the  Lennard-Jones  collision  diameter  and 
n(2-2)*  is  a  collision  integral. 

Given  the  single  component  viscosities,  we  evaluate  the  mixture  viscosity  using  the 
averaging  formula  recommended  by  Mathur  et.  al.  (5).  We  have 


(4.32) 
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Thermal  Conductivity 

The  single  component  thermal  conductivities  can  be  written  in  the  form  [6] 


y^~(ftrana^'v,trant  d"  frot^v,rot  "b  fv\b^v,vib) * 

(4.33) 

where  the  quantities  ftranatfrot  and  /OTj  are  given  by 

5/2  CVirot  A\ 
ftranB  2  \ 

(4.34) 

/,-  = fDkt  (i+-£). 

(4.35) 

f  pQkk 

foib  -  * 

(4.36) 

A=--  p^kk 

2  p* 

(4.37) 

B  =  Zro<+2  (\CT‘  +  PDkk)- 

n\Z  R  nk  ) 

(4.38) 

The  specific  heats  for  the  case  of  a  linear  molecule  take  the  familiar  forms.  We  have 


and 


Cv.rot  2 

Cvttrana  3 

(4.39) 

Cv,rot  _  , 

R 

(4.40) 
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i,vi6  —  Cv  ^R. 

(4.41) 

The  rotational  relaxation  collision  number  is  available  from  Parker  [7]  and  Brau  and 
Jonkman  [8], 

The  mixture  conductivity  is  evaluated  from  the  averaging  formula  of  Mathur  et.  al. 
(5j.  We  have 


A 


\k-i 


zLi  Xk/*k 


(4.42) 


Kinetics  Model 

The  burning  of  hydrocarbon  fuels  proceeds  through  a  complex  sequence  of  elemen¬ 
tary  reversible  reactions.  The  intermediate  species  in  such  mechanisms  are  of  extreme 
importance.  They  are  involved  in  the  chain  initiation,  chain  propagation,  chain  inhibition 
and  chain  breaking  reactions  that  describe  the  oxidation  of  a  given  fuel.  We  model  the 
kinetic  aspects  of  axisymmetric  diffusion  flames  by  allowing  the  chemical  production  and 
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destruction  of  a  species  to  occur  through  a  sequence  of  m  reversible  elementary  chemical 
reactions.  We  can  represent  these  reactions  in  the  form 


K  K 

Y*vfkZ k  **  YiujkZk,  3  =  (4.43) 

*=1  k=\ 

where  18  t^ie  stoichiometric  coefficient  of  species  k  appearing  as  a  reactant  (prod¬ 

uct)  in  reversible  reaction  j  =  1,2,  ...,m  and  where  Z ^  represents  the  chemical  symbol 
for  the  species. 

The  production  rate  w ^  for  the  species  can  be  written  in  the  form 


m 


i= l 


(4.44) 


The  function  k / (k^)  is  the  rate  constant  for  the  forward  (reverse)  path  of  reaction  j.  We 
assume  fcy  has  the  following  modified  Arrhenius  temperature  dependence 


k\  =  exp{-Ej  /  RT),  (4.45) 

with  a  similar  expression  for  k*.  The  reverse  rate  constant  can  be  written  in  terms  of  the 
forward  rate  constant  and  the  equilibrium  constants  K y  by 

kj  =  kfj/Kj‘  (4.46) 

The  pre-exponential  factor  A^,  the  temperature  exponent  fif  and  the  activation  energy  Ef 

can  be  compiled  from  published  experimental  work.  A  typical  mechanism  for  the  burning 
of  methane  in  air  is  included  in  Appendix  I. 

4.1.5.  Flame  Sheet  Starting  Estimates 

The  governing  equations  in  (4.1)-(4.12)  or  (4.16)-(4.25)  are  highly  nonlinear  and 
require  a  starting  estimate  for  the  finite  difference  method  described  in  the  next  section. 
The  determination  of  a  sufficiently  “good”  initial  solution  estimate  in  two-dimensional 
problems  can  be  difficult.  The  difficulty  is  due  to  the  exponential  dependence  of  the 
chemistry  terms  on  the  temperature  and  to  the  multidimensional  nonlinear  coupling  of 
the  fluid  and  the  thermochemistry  solution  fields.  In  previous  work,  we  focused  our  efforts 
on  the  solution  of  adiabatic  and  nonadiabatic  premixed  laminar  flames  by  adaptive  finite 
difference  methods  (9,  10,  11,  12).  In  these  problems  the  governing  differential  equations 
were  discretized  and  the  resulting  nonlinear  difference  equations  were  solved  by  Newton’s 
method.  Cubic  polynomials  and  Gaussian  shaped  profiles  were  used  as  starting  estimates 
for  the  major  and  minor  species  on  an  initial  coarse  grid.  In  adiabatic  and  nonadiabatic 
premixed  laminar  flame  problems  the  conservation  of  mass  and  momentum  reduced  to 
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the  specification  of  a  constant  mass  flow  rate  and  a  constant  thermodynamic  pressure. 
Hence,  thermochemical  considerations  played  a  more  important  role  in  these  problems 
than  did  fluid  dynamical  aspects.  This  is  rot  the  case  in  axisymmetric  laminar  diffusion 
flames.  There  is  a  strong  coupling  between  the  fluid  dynamic  and  the  thermochemistry 
solution  fields  in  these  flames.  We  have  found  that  the  solution  procedure  used  in  premixed 
laminar  flame  problems  does  not  provide  a  sufficiently  robust  or  efficient  starting  estimate 
from  which  Newton’s  method  will  converge.  In  addition,  the  relaxation  of  these  starting 
estimates  by  time-dependent  methods  to  steady-state  (or  at  least  until  the  solution  is 
within  the  convergence  domain  of  Newton’s  method)  is  extremely  slow.  The  importance  of 
these  flames  in  combustion  modeling,  however,  necessitates  the  development  of  an  efficient 
starting  procedure. 

The  burning  rate  in  a  diffusion  flame  is  controlled  by  the  rate  at  which  the  fuel  and  the 
oxidizer  are  brought  together  in  the  proper  proportions.  This  is  in  distinction  to  premixed 
flames  where  the  burning  rate  is  controlled  by  chemical  reactions.  In  the  limit  of  infinitely 
fast  kinetics,  the  fuel  and  the  oxidizer  are  separated  by  a  thin  exothermic  reaction  zone. 
In  this  zone  the  fuel  and  the  oxidizer  are  in  stoichiometric  proportion  and  the  temperature 
and  products  of  combustion  are  maximized.  In  such  an  ideal  situation,  no  oxidizer  is 
present  on  the  fuel  side  and  no  fuel  is  present  on  the  oxidizer  side.  The  fuel  and  oxidizer 
diffuse  towards  the  reaction  zone  as  a  result  of  concentration  gradients  in  the  flow.  In 
diffusion  flames  of  practical  interest,  the  oxidation  of  the  fuel  to  form  intermediates  and 
products  proceeds  through  a  detailed  kinetics  mechanism.  In  these  problems  combustion 
takes  place  at  a  finite  rate  and  some  fuel  and  oxidizer  co-exist  on  either  side  of  the  reaction 
zone.  Nevertheless,  the  use  of  a  thin,  infinitely  fast,  global  reaction  model  is  a  natural 
starting  point  for  the  determination  of  a  “good”  initial  solution  estimate  for  our  finite  rate 
axisymmetric  model  (see  also  Burke  and  Schumann  [13],  Mitchell  [2],  Mitchell  et  al.  [14], 
Smooke  et  al.  [15],  and  Keyes  and  Smooke  [16]). 

We  assume  that  the  fuel  and  the  oxidizer  obey  a  single  overall  irreversible  reaction  of 
the  type 

Fuel  ( F )  -f-  Oxidizer  (X)  — »  Products  (P),  (4.47) 

in  the  presence  of  an  inert  gas  (N).  We  have 

upF  +  f  t'pP,  (4.48) 

where  up,ux  and  up  are  the  stoichiometric  coefficients  of  the  fuel,  the  oxidizer  and  the 
product,  respectively.  In  addition,  we  neglect  thermal  diffusion  and  assume  that  cp  and 
cpk  are  constant  and  that  the  ordinary  mass  diffusion  velocities  can  be  written  in  terms  of 
Fick’s  law.  With  these  approximations  we  can  write 

Stream  Function: 


Vorticity: 


(4.49) 
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-h  K  (’"))  -  i  K  (’•))  >-'4 


Species: 


and 

Energy: 


where 


~lL  (rpDf'lff)  ~  rWFuF^  =  0. 

-[£(■>£)*£(*£)]- 1  (**£) 

(jpDx~lhf)  ~  rWXuX^  =  0. 

-[£(>>SK(^)] -£(•'<?) 

-£  ('^p£r)  ‘ '“w»  - 

-[£«)*£(*--S)]-£K^) 

)-£(•>!) 

Wfi/fhf  -f  W) ~  Wpvphp 


— r- 


•u>  =  0, 


u;  =  -^L  = 


up  t'x  vp 

is  the  rate  of  progress  of  the  reaction. 

If  we  introduce  the  heat  release  per  unit  mass  of  the  fuel  Q  where 

^  l  ,  WXUX,.  WPvpL 
Q  =  hp+  r  hX  ~  777 —rhP* 


Wpvp 

and  if  we  assume  that  the  Lewis  numbers 

A 


Wpup 


Lejr  = 


pDpCp 


,Lex  = 


pDxCp  ’ 


0,  (4.50) 

(4.51) 

(4.52) 

(4.53) 

(4.54) 

(4.55) 

(4.56) 

(4.57) 

(4.58)o 
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T.  A  T-  _  A 

Lcp  —  n  »  Lc/y  —  n  y 

pDpcr  pDtfCp 

(4.59)6 

are  all  equal  to  one,  then  each  of  the  Shvab-Zeldovich  variables 

ZF  =  YF-  YFq  +  |(T  -  To), 

(4.60) 

Zx  =  yx~yXo  +  ^wwxfIxf(t-to), 

(4.61) 

Zr=YP-YPo-'^JF(T-T0), 

(4.62) 

ZN  =  Yn  ~  YNq, 

(4.63) 

satisfies  the  differential  equation 


(4.64) 


One  can  show  that  all  of  the  Z *  are  proportional  to  each  other  and  to  a  conserved  scalar 
5  which  satisfies  an  equation  similar  in  form  to  (4.64). 

To  complete  the  specification  of  the  starting  estimate,  we  must  be  able  to  recover 
the  temperature  and  the  major  species  profiles  from  the  conserved  scalar.  Of  critical 
importance  to  this  procedure  is  an  estimate  of  the  location  of  the  flame  front.  In  the 
Shvab-Zeldovich  formulation  fuel  and  oxidizer  cannot  co-exist.  Hence,  on  the  fuel  side  of 
the  flame  Y\  —  0  and  on  the  oxidizer  side  Yp  —  0.  If  we  denote  variables  at  the  flame 
front  with  the  subscript  /,  then  it  can  be  shown  (see  also  Keyes  and  Smooke  [16])  that, 
for  a  fixed  value  of  the  axial  coordinate  z,  the  location  of  the  flame  front  is  defined  such 
that 

■S(r/)l/.I<J  .  =  S,  =  (4.65) 

The  location  of  the  flame  front  can  be  obtained  by  solving  (4.65)  at  each  axial  coordinate 
level. 

If  we  utilize  the  proportionality  of  the  Z^’a  to  S  along  with  the  expressions  in  (4.49)- 
(4.63),  we  can  derive  expressions  for  the  temperature  and  species  on  the  fuel  and  oxidizer 
sides  of  the  flame.  On  the  fuel  side,  we  have 


T  =  r/s+[r0  +  yXo9^](i-5), 

(4.66) 

rF-rns  +  YJbZ£'(s  i), 

(4.67) 
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and 


YX  =  0, 


YP  =  YXo 


«W(1 


-5), 


(4.68) 

(4.69) 


=  ^OVoi1  ~  5) + YNfS. 

On  the  oxidizer  side,  we  have 


T  =  T0{1  -S)  + 


i?r*  +  r' 


s. 


(4.70) 


(4.71) 


Yf  =  0, 

^  =  ^(1-5)-^,^. 


fVpt/p 

Kp  -  Ww  F'S' 


(4.72) 

(4.73) 

(4.74) 


and 

Ytl  =  YKo{1  -  s)  +  Yn,S-  MS) 

The  equations  in  (4.49),  (4.50)  and  (4.64)  are  solved  for  the  stream  function,  the 
vorticity,  the  temperature  and  the  major  species  in  the  flame.  For  a  given  profile  of  the 
conserved  scalar,  we  solve  (4.65)  for  the  location  of  the  flame  front  on  each  axial  level. 
We  then  utilize  the  relations  in  (4.66)-(4.75)  to  obtain  expressions  for  T,Yf,  Yq,  Yp  and 
Yff.  The  recovered  temperature  profile  is  used  in  the  ideal  gas  law  to  evaluate  the  density. 
The  temperature  is  also  needed  in  forming  the  viscosity  and  the  diffusion  coefficient  [D^ 
is  replaced  by  D).  If  we  introduce  the  Prandtl  number 


Pr  = 


HCP 

A  ’ 


(4.76) 


and  recall  that  all  of  the  Lewis  numbers  are  equal  to  one,  we  can  write 


(4.77) 


where  (Pr)rey  is  a  reference  Prandtl  number.  Specifically,  we  use  an  approximate  value 
for  air,  (Pr)r<^  =  0.75.  Hence,  determination  of  pD  is  reduced  to  the  specification  of  a 
transport  relation  for  the  viscosity.  We  use  the  simple  power  law 


(4.78) 


where  r  -  0.7,  T0  =  298  K  and  f. i0  =  1.85  x  10~4  gm/cm-sec  is  again  a  reference  value 
for  air  (17).  The  temperature  exponent  was  determined  by  fitting  the  equation  in  (4.78) 
to  the  mixture  viscosity  and  temperature  data  of  a  representative  one-dimensional  finite 
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rate  chemistry  calculation.  The  scaled  heat  release  parameter  Q/cp  can  be  determined 
from  an  estimate  of  the  peak  temperature  ( e.g.,  from  an  experiment)  or  from  the  heat  of 
combustion  of  the  system  under  consideration  and  a  representative  heat  capacity. 

4.1.4.  Computational  Approach 

Most  of  the  detailed  chemistry,  computational  combustion  studies  that  have  appeared 
in  the  literature  have  focused  on  essentially  one-dimensional  configurations,  i.e.,  freely 
propagating  or  burner-stabilized  premixed  flames  and  counterflow  premixed  or  diffusion 
flames  (Spalding  (18),  Adams  and  Cook  [19],  Dixon-Lewis  (20,  21),  Spalding  et  al.[22j, 
Wilde  (23),  Bledjian  (24),  Margolis  (25),  Wamatz  (26,  27,  28),  Westbrook  and  Dryer  (29, 
30],  Coffee  and  Heimerl  (31,  32),  Miller  et  al.  (9,  10,  33),  Smooke  (llj,  Smooke  et  al.  (34, 
12,  35,  36),  Hahn  and  Wendt  (37),  Sato  and  Tsuji  (38),  Dixon-Lewis  et  al.  (39),  Giovangigli 
and  Smooke  (40)).  The  interaction  of  heat  and  mass  transfer  and  chemical  reaction  in 
practical  combustion  systems,  however,  requires  a  multidimensional  study. 

Our  goal  is  to  obtain  a  discrete  solution  of  the  governing  equations  in  two  dimensions 
on  the  mesh  M2  the  initial  nodes  of  which  are  formed  by  the  intersection  of  the  lines  of 
the  mesh  Mr 

Mr  =  (0  =  r0  <  rj  <  . . .  <  r, . . .  <  rwr  =  R0],  (4-79) 

and  the  mesh  Mx 

M,  =  {0  =  zq  <  z\  <  ...  <  Zj  ...  <  zmm  =  Z }.  (4.80) 

Computationally,  we  combine  a  steady-state  and  a  time-dependent  solution  method.  A 
time-dependent  approach  is  used  to  help  obtain  a  converged  numerical  solution  on  an 
initial  coarse  grid  using  the  flame  sheet  starting  estimate.  Grid  points  are  then  inserted 
adaptively  and  the  steady-state  solution  procedure  is  used  to  complete  the  problem. 

Newton’s  Method 

We  approximate  the  spatial  operators  in  the  governing  partial  differential  equations  by 
finite  difference  expressions.  Diffusion  terms  are  approximated  by  centered  differences  and 
convective  terms  by  upwind  approximations.  The  problem  of  finding  an  analytic  solution 
of  the  equations  is  then  converted  into  one  of  finding  an  approximation  to  this  solution  at 
each  point  (r^Zj)  of  the  mesh  in  two  dimensions.  With  the  difference  equations  written 
in  residual  form,  we  seek  the  solution  U *  of  the  system  of  nonlinear  equations 

F(U)  =  0.  (4.81) 

For  an  initial  solution  estimate  U°  which  is  sufficiently  close  to  U*,  the  system  of  nonlinear 
equations  in  (4.81)  can  be  solved  by  Newton’s  method.  This  leads  to  the  iteration 

J(Un)(C/n+1  -  Un)  =  -XnF[Un)t  n  =  0,1,2,....  (4.82) 

J(Un)  =  dF(Un)/dU  is  the  Jacobian  matrix  and  Xn  (0  <  A  <  1)  is  the  damping 
parameter  (41). 

We  point  out  that  with  the  spatial  discretizations  used  in  forming  (4.81),  the  Jacobian 
matrix  in  (4.82)  can  be  written  in  block  nine  diagonal  form.  For  problems  involving 


15 


detailed  transport  and  complex  chemistry,  it  is  often  more  efficient  to  evaluate  the  Jacobian 
matrix  numerically  as  opposed  to  analytically.  The  numerical  procedure  we  implement 
extends  the  ideas  outlined  by  Curtis,  Powell  and  Reid  (42).  We  form  several  columns  of 
the  Jacobian  simultaneously  using  vector  function  evaluations  and  the  Jacobian's  given 
sparsity  structure.  If  to  each  column  of  the  Jacobian  we  associate  the  *  and  j  values  of 
the  node  corresponding  to  the  column’s  diagonal  block,  then  all  columns  of  the  Jacobian 
having  the  same  value  of  the  parameter 


a  =  (»  4-  3j)mod  9, 


(4.83) 


can  be  evaluated  simultaneously.  Ideas  along  these  lines  have  also  been  explored  by 
Newsam  and  Ramsdell  [43]  and  Coleman  and  More  [44].  Once  the  Jacobian  is  formed 
we  solve  the  Newton  equations  with  a  block-line  SOR  method.  The  Newton  iteration 
continues  until  the  size  of  ||£/n+1  —  Un ||2  is  reduced  appropriately. 

Adaptive  Gridding 

The  solutions  of  the  governing  equations  in  the  axisymmetric  problems  contain  regions 
in  each  coordinate  direction  in  which  the  dependent  variables  exhibit  high  spatial  activity 
(steep  fronts  and  sharp  peaks).  Efficient  solution  of  these  problems  requires  that  the  high 
activity  regions  be  resolved  adaptively.  Techniques  that  attempt  to  equidistribute  positive 
weight  functions  have  been  used  with  a  great  deal  of  success  in  premixed  and  counterflow 
flame  problems.  Flames  with  30  to  40  chemical  species  and  over  100  chemical  reactions 
can  be  solved  efficiently  by  adaptively  placing  grid  points  in  the  high  activity  regions.  We 
follow  a  similar  approach  in  the  axisymmetric  problem. 

Adaptive  mesh  refinement  in  two  dimensions  can  proceed  along  several  different  paths. 
The  simplest  procedure  involves  determining  the  grid  points  of  the  mesh  M2  by  equidis- 
tributing  positive  weight  functions  over  mesh  intervals  in  both  the  r  and  z  directions  (see, 
e.g.,  Kautsky  and  Nichols  [45]  and  Russell  [46]).  Specifically,  we  attempt  to  equidistribute 
the  mesh  Mr  with  respect  to  the  non-negative  function  Wr  and  constant  Cr  for  each  of  the 
Ms  +  1  horizontal  grid  lines.  We  write 


/; 


«+i 


Wr  dr  <  Cr,  *  =  0,  l,...,Afr  —  1, 


(4.84) 


for  j  =  0, 1, . . . ,  Similarly,  we  attempt  to  equidistribute  the  mesh  M*  with  respect  to 
the  non-negative  function  Ws  and  constant  Cg  for  each  of  the  Mr  +  1  vertical  grid  lines. 
We  have 

“*j+i 

Wa  dz<Ca ,  j  =  0,1,..., Af*  —  1,  (4.85) 


f 

J 


for  «  =  0, 1,. . . , Mr. 

In  implementing  the  two-dimensional  adaptive  grid  strategy,  we  first  solve  the  bound¬ 
ary  value  problem  on  an  initial  coarse  grid.  We  then  test  the  inequality  in  (4.84)  one  r 
subinterval  at  a  time  for  all  the  j  grid  lines  and  all  the  dependent  solution  components.  If 
the  inequality  is  not  satisfied,  a  grid  point  is  inserted  at  the  midpoint  of  the  r  subinterval  in 
question  for  j  =0,1,...,  Mg.  Once  this  procedure  has  been  carried  out  in  the  r  direction, 
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we  reverse  the  process  and  begin  again  in  the  z  direction.  This  procedure  produces  an 
orthogonal  tensor  product  grid  in  which  the  coordinate  lines  connect  opposite  boundaries 
of  the  computational  domain.  The  weight  functions  in  the  equidistribution  procedure  are 
chosen  such  that  the  grid  points  are  placed  in  regions  of  high  spatial  activity  with  the 
goal  of  reducing  the  local  discretization  error.  We  use  a  combination  of  first  and  second 
derivatives  of  the  solution  profiles  (see  also  Smooke  [11,  47),  Giovangigli  and  Smooke  [40)). 
The  particular  combinations  of  function  and  slope  and  the  values  of  Cr  and  Ca  can  be 
changed  to  produce  a  solution  to  a  desired  level  of  accuracy. 

Coarse  to  Fine  Grid  Methodology 

The  formation  of  the  Jacobian  and  its  partial  factorization  in  the  block-line  SOR 
method  accounts  for  a  substantial  part  of  the  cost  of  the  diffusion  flame  calculation.  As 
a  result,  the  use  of  a  modified  Newton  method  in  which  the  Jacobian  is  re-evaluated 
periodically  is  indicated.  The  immediate  implication  of  applying  the  modified  Newton 
method  is  that  the  partial  factorization  of  the  Jacobian  can  be  stored  and  each  modified 
Newton  iteration  can  be  obtained  by  performing  relatively  inexpensive  block-line  SOR 
back  substitutions.  The  problem  one  faces  when  applying  the  modified  method  is  how  to 
determine  whether  the  rate  of  convergence  is  fast  enough.  If  the  rate  is  too  slow  we  want 
to  change  back  to  a  full  Newton  method  and  make  use  of  new  Jacobian  information.  If 
the  rate  of  convergence  is  acceptable,  we  want  to  continue  performing  modified  Newton 
iterations. 

In  our  adaptive  grid  strategy,  the  equidistribution  condition  is  checked  one  mesh 
interval  at  a  time  and  grid  points  are  added  appropriately.  The  coarse  grid  solution  is  then 
interpolated  linearly  onto  the  new  finer  grid.  The  interpolated  result  serves  as  an  initial 
solution  estimate  for  the  iteration  procedure  on  the  finer  grid.  The  process  is  continued 
on  successively  finer  and  finer  grids  until  several  termination  criteria  are  satisfied.  We 
anticipate  that  as  the  size  of  the  mesh  spacing  gets  smaller,  the  interpolated  solution 
should  become  a  better  starting  estimate  for  Newton’s  method  on  the  next  finer  grid.  For 
a  class  of  nonlinear  boundary  value  problems,  Smooke  and  Mattheij  [48)  have  shown  that 
there  exists  a  critical  mesh  spacing  such  that  the  interpolated  solution  lies  in  the  domain 
of  convergence  of  Newton’s  method  on  the  next  grid.  As  a  result,  the  hypotheses  of  the 
Kantorovich  theorem  [49)  are  satisfied  and  the  sequence  of  successive  modified  Newton 
iterates  can  be  shown  to  satisfy  a  recurrence  relation  scaled  by  the  first  Newton  step 
[50].  As  a  result,  if  in  the  course  of  a  calculation,  we  determine  that  the  size  of  the 
(n-f- 1)8^  modified  Newton  step  is  larger  than  the  value  predicted  by  the  theorem,  we  form 
a  new  Jacobian  and  restart  the  iteration  count. 

Time-Dependent  Starting  Estimates 

The  coarse  to  fine  grid  strategy  and  the  flame  sheet  starting  estimate  helps  elimi¬ 
nate  many  of  the  convergence  difficulties  associated  with  solving  the  governing  equations 
directly.  Nevertheless,  to  obtain  a  starting  estimate  on  the  initial  grid  that  lies  in  the 
convergence  domain  of  Newton’s  method  we  apply  a  time-dependent  iteration  to  the  flame 
sheet  solution.  We  remark  that  fundamentally  there  are  two  mathematical  approaches  for 
solving  flame  problems  -  one  uses  a  transient  method  and  the  other  solves  the  steady-state 
boundary  value  problems  directly.  Generally  speaking,  the  transient  methods  are  robust 
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but  computationally  inefficient  compared  to  the  boundary  value  methods,  which  are  effi¬ 
cient  but  have  less  desirable  convergence  properties.  Most  of  the  numerical  techniques  that 
have  been  used  to  solve  one-dimensional  flame  problems  have  employed  a  time-dependent 
method.  Variations  of  this  approach  have  been  considered  by  a  variety  of  researchers  (see, 
e.g.,  Spalding  (18],  Adams  and  Cook  (19],  Dixon-Lewis  (20,  21],  Spalding  et  al.  [22],  Bled- 
jian  [24],  Margolis  (25],  Warnatz  (26,  27,  28],  Westbrook  and  Dryer  (29,  30],  Coffee  and 
Heimerl  (31,  32]).  In  these  methods,  the  original  nonlinear  two-point  boundary  value  prob¬ 
lem  is  converted  into  a  nonlinear  parabolic  mixed  initial-boundary  value  problem.  This 
is  accomplished  by  appending  the  term  d[*)/dt  to  the  left-hand  side  of  the  conservation 
equations.  This  same  procedure  can  be  employed  in  our  two-dimensional  calculations.  We 
obtain 

^  =  F(U),  (4.86) 

with  appropriate  initial  conditions.  If  the  time  derivative  is  replaced,  for  example,  by  a 
backward  Euler  approximation,  the  governing  equations  can  be  written  in  the  form 

/(t/n+1)  =  F{Un+l)  -  ~  =  0,  (4.87) 

where  for  a  function  g(t)  we  define  gn  =  g{tn)  and  where  the  time  step  rn+1  =  tn+1  —  tn. 
At  each  time  level  we  must  solve  a  system  of  nonlinear  equations  that  look  very  similar 
to  the  nonlinear  equations  in  (4.81).  Newton’s  method  can  again  be  used  to  solve  this 
system.  The  important  difference  between  the  system  in  (4.81)  and  (4.87)  is  that  the 
diagonal  of  the  steady-state  Jacobian  is  weighted  by  the  quantity  l/rn+1.  This  produces  a 
better  conditioned  system  and  the  solution  from  the  time  step  ordinarily  provides  an 
excellent  starting  guess  to  the  solution  at  the  (n  +  l)8*  time  level.  The  work  per  time  step 
is  similar  to  that  for  the  modified  Newton  iteration,  but  the  timelike  continuation  of  the 
numerical  solution  produces  an  iteration  strategy  that  will,  in  general,  be  less  sensitive  to 
the  initial  starting  estimate  than  if  Newton’s  method  were  applied  to  (4.81)  directly.  As 
a  result,  when  we  ultimately  implement  Newton’s  method  on  the  steady-state  equations 
directly,  we  obtain  a  converged  numerical  solution  with  only  a  few  additional  iterations. 
This  time-dependent  starting  procedure  can  also  be  used  on  grids  other  than  the  initial 
one.  The  size  of  the  time  steps  are  chosen  by  monitoring  the  local  truncation  error  of  the 
time  discretization  process  (see  also,  Smooke  et  al.  [35]). 

4.1.5.  Numerical  Results  —  Laminar  Flames 

In  this  section  we  apply  the  flame  sheet  starting  estimate  and  the  computational 
method  discussed  in  the  previous  section  to  both  a  confined  and  an  unconfined  methane- 
air,  axi8ymmetric  laminar  diffusion  flame.  Detailed  transport  coefficients  and  a  42  reaction, 
fifteen  species  reaction  mechanism  (see  Table  I,  Miller  et  al.  [33]  and  Smooke  et  al.  [36] 
were  used  in  the  calculations. 

Confined  Flame 

The  first  flame  we  consider  is  a  confined  methane-air  diffusion  flame  studied  exper¬ 
imentally  by  Mitchell  [2].  The  experimental  configuration  is  such  that  the  radius  of  the 
inner  fuel  jet  Rj  =  0.635  cm,  the  radius  of  the  outer  oxidizer  jet  Rq  =  2.54  cm  and  the 
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length  of  the  tubular  pyrex  shield  Z  =  30.0  cm.  Fuel  is  introduced  through  the  center 
tube  and  air  through  the  outer  coflow.  The  boundary  conditions  at  the  inlet  are  given  by 

Inlet  ( z  =  0) : 

imrr - 


T  =  298  K, 

YCH4  =  1.0,  Yk  =  0,  k  *  CHa>  (4.88) 

vr  =  0.0  cm/sec, 
v„  =  4.5  cm/sec, 

Rj  <  r  <  Rq 


T  =  298  K, 

Y0j  =  .232,  YNi  =  .768,  Yk  =  0,  k  *  02tN2t  (4.89) 

vr  =  0.0  cm/sec, 
vs  =  9.88  cm/sec. 

The  appropriate  boundary  conditions  for  the  stream  function-vorticity  system  can  be 
formed  using  the  definition  in  (4.13)  along  with  (4.22)-(4.23).  The  shield  temperature 
is  kept  constant  at  298  K. 

The  flame  sheet  model  provided  initial  solution  profiles  for  the  stream  function,  the 
vorticity,  the  temperature  and  the  major  species,  i.e.,  CH4,  02,  N2,  C02  and  H20.  The 
starting  estimates  for  the  minor  species  in  the  full  chemistry  solution  were  approximated 
by  Gaussian  profiles  that  were  centered  at  the  location  of  the  flame  sheet  on  each  axial 
level.  They  had  peak  heights  of  at  most  a  few  percent.  To  conserve  mass  in  the  starting 
estimate,  the  N2  mass  fraction  was  reduced  accordingly.  The  flame  sheet  starting  estimate 
required  approximately  150  adaptive  time  steps  and  five  Newton  iterations  to  converge. 
Once  the  flame  sheet  estimate  was  calculated,  we  solved  the  full  set  of  governing  equations 
in  a  two-step  procedure.  We  first  determined  a  solution  to  the  stream  function,  vorticity 
and  species  equations  (4.16)— (4.18)  based  on  the  flame  sheet  temperature  profile.  This 
fixed  flame  sheet  temperature  solution  {Tout)  was  then  used  as  input  to  the  full  fluid 
dynamic-thermochemistry  model  (4.16)-(4.19)  in  which  the  energy  equation  was  included 
{T[ff).  This  procedure  helped  to  reduce  both  convergence  difficulties  and  the  total  CPU 
time  and  is  similar  to  the  two-pass  solution  method  used  in  the  solution  of  adiabatic 
premixed  laminar  flames  (12]  and  counterflow  diffusion  flames  [16]. 

The  flame  sheet  and  the  first  Tqut  calculation  were  performed  on  a  40  x  28  tensor 
product  grid.  One  hundred  fifty-five  adaptive  time  steps  were  required  to  reduce  the 
norm  of  the  Tout  steady-state  residuals  below  1.0  x  10~*.  This  was  sufficient  to  bring 
the  numerical  solution  within  the  convergence  domain  of  Newton’s  method.  After  the 
time  steps,  Newton’s  method  converged  with  only  seven  iterations.  Once  this  solution 
was  obtained,  the  mesh  was  refined  and  a  solution  was  calculated  on  a  finer  grid.  This 
procedure  was  continued  until  a  refined  75  x  41  grid  was  obtained.  The  refined  fixed 
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temperature  solution  was  then  used  as  the  starting  estimate  for  the  complete  fluid  dynamic- 
thermochemistry  solution.  As  the  computational  mesh  was  refined,  Newton’s  method 
typically  converged  with  a  smaller  number  of  time  steps  than  on  the  coarser  grids.  The 
mesh  spacing  was  such  that  32000  equispaced  points  would  have  been  needed  to  obtain 
comparable  accuracy.  The  total  CPU  time  for  the  entire  procedure  was  approximately 
150  hours  on  an  FPS-264  which  is  consistent  with  a  single  one-dimensional  counterflow 
diffusion  flame  calculation  that  requires  between  2-3  hours  on  the  same  machine. 

Mitchell  [2]  has  measured  radial  temperature  profiles  at  several  heights  above  the 
burner  using  three  mil  Pt  vs.  Pt-13%  Rh  thermocouples.  The  measurements  were  cor¬ 
rected  for  radiative  losses.  Major  species  were  sampled  with  quartz  microprobes  and  then 
analyzed  with  a  gas  chromatograph.  The  probe  had  a  maximum  outside  diameter  of  0.75 
mm  and  was  tapered  to  0.15  mm  at  the  tip.  The  inner  diameter  was  approximately  50/x. 
In  Figure  3  we  compare  the  radial  experimental  and  calculated  temperature  profiles  at 
a  height  of  1.2  cm  above  the  burner.  We  observe  excellent  agreement  between  the  two 
profiles  from  the  axis  of  symmetry  to  the  low  temperature  coflow  region.  In  Figures  4  and 
5  we  compare  radial  experimental  and  computational  profiles  for  the  major  species  in  the 
flame  ( CH4t02,N2,H20,CC>2,C0  and  Jfy)  at  a  height  of  1.2  cm  above  the  burner.  The 
agreement  among  all  the  computational  and  experimental  profiles  is  very  good.  Several 
features  in  Figures  4  and  5  are  worth  noting,  however.  From  Figure  4  it  is  clear  that 
the  methane  diffuses  to  the  reaction  zone  where  it  is  completely  consumed.  The  oxygen 
concentration  outside  the  flame  region  is  near  its  inlet  value  and  then  drops  nearly  to  zero 
in  the  reaction  zone.  It  then  slightly  increases  as  the  symmetric  axis  is  approached.  This 
type  of  oxygen  profile  results  near  the  flame  base  because  the  temperatures  are  low  enough 
to  allow  some  oxygen  to  penetrate  the  flame.  In  the  radial  direction  water  and  carbon 
dioxide  maximize  in  the  region  of  the  peak  temperature.  In  the  lower  regions  of  the  flame 
the  carbon  monoxide  and  hydrogen  profiles  first  increase  with  distance  from  the  symmetric 
axis  then  proceed  through  a  maximum  and  finally  decrease  to  zero  in  the  reaction  zone. 
Similar  profiles  for  the  temperature  and  major  species  are  illustrated  in  Figures  6-8  at  a 
height  of  2.4  cm  above  the  burner.  We  note  that  the  experimental  temperature  profile  is 
somewhat  narrower  than  the  computed  one  but  the  peak  values  and  the  trend  in  the  major 
species  profiles  is  again  excellent.  Experimentally,  temperature  profiles  were  not  obtained 
from  2.5  cm  to  5.0  cm  above  the  burner  due  to  the  tendency  of  carbon  particles  to  collect 
on  the  thermocouple  bead. 

A  more  global  picture  of  this  flame  can  be  obtained  by  plotting  the  temperature  and 
species  contours  versus  the  independent  spatial  coordinates.  In  Figures  9-11  we  illustrate 
the  temperature  isotherms  and  the  methane  and  oxygen  isopleths,  respectively,  as  functions 
of  the  axial  and  radial  coordinates.  We  note  immediately  the  high  temperature  region 
extending  from  the  boundary  of  the  fuel  and  oxidizer  jets  to  the  axis  of  symmetry.  The 
figure  also  points  out  the  extremely  high  temperature  gradients  directly  above  the  burner 
inlet.  The  temperature  rises  from  298  K  to  nearly  2000  K  in  approximately  0.8  mm.  It 
is  in  this  region  that  the  fuel  and  the  oxidizer  first  meet  in  stoichiometric  proportion.  In 
particular,  it  is  clear  from  these  figures  that  combustion  occurs  only  in  a  thin  region  above 
the  inlet.  The  methane  and  oxygen  co-exist  in  only  a  very  small  region.  Most  of  the 
methane  disappears  within  1.0  cm  of  the  fuel  jet.  The  resulting  heat  release  produces  an 
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extremely  rapid  rise  in  the  temperature.  As  one  moves  in  the  direction  of  increasing  r,  we 
observe  that  the  temperature  rapidly  decreases  and  ultimately  approaches  its  inlet  value. 
In  addition,  if  we  define  the  flame  height  as  the  first  social  location  where  the  maximum 
temperature  occurs  on  the  axis  of  symmetry,  we  obtain  a  flame  height  of  approximately 
7.0  cm. 

From  a  structural  viewpoint,  methane  is  the  simplest  hydrocarbon-it  is  the  only  hy¬ 
drocarbon  without  a  carbon-carbon  bond.  As  a  result,  its  oxidation  differs  significantly 
from  other  hydrocarbon  fuels.  Westbrook  and  Dryer  [51]  and  Warnatz  [52}  have  postu¬ 
lated  that  the  oxidation  of  methane  occurs  through  roughly  two  parallel  paths.  In  one 
path  carbon-hydrogen  bonds  are  broken  to  form  methyl  radicals.  The  methyl  radicals  are 
then  oxidized  to  methoxy  radicals  and/or  formaldehyde.  This  is  followed  by  the  formation 
of  the  formal  group  which  is  then  oxidized  to  form  carbon  monoxide.  In  the  second  path, 
methyl  radical  recombination  is  followed  by  the  oxidation  of  the  resulting  Ci  species.  In 
our  calculations,  the  former  path  was  chosen  and  the  latter  path  was  neglected. 

In  Figures  12-15  we  illustrate  isopleths  for  ^0,00,  Hi  and  CO 2,  respectively.  Sim¬ 
ilar  plots  for  CH$tOH  and  CH2O  are  illustrated  in  Figures  16-18.  It  is  clear  from  these 
figures  that  large  quantities  of  H^O^CO  and  H2  are  produced  soon  after  the  methane 
has  been  consumed.  It  is  in  this  region  that  the  methane  is  attacked  by  0,H  and  OH 
radicals  and  CH3  is  formed.  In  this  region  only  small  amounts  of  OH ,  H  and  O  exist  (see, 
e.g .,  Figure  17)  due  to  the  high  affinity  of  methane  for  these  radicals.  The  peak  values 
of  both  CH$  and  CH2O  also  occur  after  the  methane  has  disappeared.  The  oxidation  of 
CH2O  to  HCO  and  the  subsequent  formation  of  CO  occurs  in  regions  of  high  methyl  and 
formaldehyde  concentrations. 

The  oxidation  of  CO  to  CO2  proceeds  primarily  via  the  reaction 

CO  +  OH  -►  C02  +  H. 

Hence,  the  rate  of  CO  oxidation  depends  on  the  availability  of  OH  radicals.  However,  as 
Westbrook  and  Dryer  point  out  [51]  the  presence  of  most  hydrocarbon  species  inhibits  the 
oxidation  of  CO.  This  can  be  attributed  to  the  fact  that  the  rate  of  the  reaction 

H  +  02->  OH  +  0, 

is  considerably  smaller  than  the  reaction  rates  of  H  atoms  with  hydrocarbon  species  and 
the  rate  of  the  CO  oxidation  reaction  is  also  smaller  than  the  reaction  rates  of  hydrocarbon 
species  with  OH.  As  a  result,  small  quantities  of  hydrocarbons  can  effectively  restrict  the 
oxidation  of  CO  to  CO2 .  Although  carbon  monoxide  and  hydrogen  are  found  during 
the  oxidation  of  the  hydrocarbon  species,  it  is  not  until  after  the  hydrocarbons  and  the 
hydrocarbon  fragments  have  been  consumed  that  the  OH  level  rises  and  CO2  is  formed. 
In  Figure  17  we  observe  that  in  the  axial  direction  the  OH  radical  pool  increases  after 
the  disappearance  of  the  methane  and  the  formation  of  CO.  The  CO  is  then  oxidized  to 
form  CO2 •  Hence,  in  Figure  15  we  observe  that  carbon  dioxide  forms  downstream  of  the 
regions  of  high  CH$,CH20  and  CO  concentration. 

From  a  fluid  dynmaic  viewpoint,  we  find  the  velocities  along  the  centerline  to  be  higher 
than  the  values  one  would  expect  due  to  the  effects  of  natural  convection.  In  particular, 
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the  stream  function  isopleths  in  Figure  19  illustrate  two  large  recirculation  cells  that  are 
established  between  the  hot  surface  of  the  flame  and  the  cooler  shield.  Air  is  entrained  into 
the  system  at  the  shield  outlet  to  balance  the  momentum  of  the  inlet  fuel  and  air  streams 
along  with  the  frictional  losses  at  the  shield  wall.  The  presence  of  these  recirculation  cells 
reduces  the  total  area  available  for  the  flow  of  the  combustion  gases  and  hence  the  velocities 
are  increased  due  to  the  combined  effects  of  natural  convection  and  a  reduced  flow  area. 

It  is  worthwhile  to  reconsider  Figure  8  in  light  of  the  recirculation  cells.  We  note  that 
the  water  vapor  profile  in  Figure  8  asymptotes  to  a  value  of  about  one  mole  percent.  (At 
a  height  of  5.0  cm  above  the  burner  this  value  asymptotes  to  about  three  mole  percent). 
This  is  due  to  the  fact  that  moist  air  was  entrained  into  the  system  at  the  shield  outlet. 
The  moisture  content  of  the  air  in  the  exit  region  of  the  shield  was  higher  than  in  other 
portions  of  the  laboratory.  As  a  result,  following  Mitchell’s  [2]  measurements,  the  inflow 
boundary  conditions  at  the  top  of  the  shield  in  the  recirculation  zone  were  adjusted  to 
account  for  the  higher  water  vapor  concentrations  in  the  entrained  air. 

Unconfined  Flame 

The  second  flame  we  consider  is  an  unconfined  methane-air  diffusion  flame.  The 
experimental  configuration  is  such  that  the  radius  of  the  inner  fuel  jet  Rj  =  0.2  cm  and 
the  radius  of  the  outer  oxidizer  jet  Rq  =  2.54  cm.  Fuel  is  introduced  through  the  center 
tube  and  air  through  the  outer  coflow.  The  boundary  conditions  at  the  inlet  are  given  by 

Inlet  ( z  —  0): 

- 


T  =  298  K, 

YCH<  =  1.0,  n=0,  M  CH4,  (4.90) 

Vf  —  0.0  cm/sec, 
v,  =  5.0  cm/sec, 

Rj  <  r  <  Rq 


T  =  298  K, 

V0j  =  .232,  YNj  =  .768,  Yk  =  0,k  ±  02,N2,  (4.91) 

vr  =  0.0  cm/sec, 
v,  =  25.0  cm/sec. 

The  appropriate  boundary  conditions  for  the  stream  function-vorticity  system  can  be 
formed  using  the  definition  in  (4.13)  along  with  (4.22)-(4.23). 

The  flame  sheet  model  again  provided  initial  solution  profiles  for  the  stream  function, 
the  vorticity,  the  temperature  and  the  major  species,  ».e.,  CH4 ,  02 ,  N2,  C02  and  H20 
and  the  starting  estimates  for  the  minor  species  in  the  full  chemistry  solution  were  also 
approximated  by  Gaussian  profiles  that  were  centered  at  the  location  of  the  flame  sheet  on 
each  axial  level.  The  two-step  solution  procedure  described  in  the  confined  flame  problem 
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was  used  in  solving  the  complete  set  of  governing  equations.  The  flame  sheet  and  the  first 
Tout  calculation  were  performed  on  a  31  x  25  tensor  product  grid.  The  Tout  solution 
was  refined  four  times  until  a  solution  was  obtained  on  a  65  x  38  grid.  The  total  CPU  time 
was  comparable  to  the  time  required  to  obtain  a  solution  for  the  confined  flame. 

From  a  chemical  viewpoint,  the  flame  is  structurally  similar  to  the  confined  flame 
discussed  previously.  In  Figure  20  we  illustrate  the  temperature  isotherms  for  this  system. 
We  again  observe  extremely  high  temperature  gradients  directly  above  the  burner  inlet. 
The  methane  and  oxygen  co-exist  in  only  a  very  small  region.  Most  of  the  methane  disap¬ 
pears  within  1.0  mm  of  the  fuel  jet.  Similarly,  water,  carbon  monoxide  and  hydrogen  are 
produced  soon  after  the  methane  has  been  consumed.  The  oxidation  of  CO  to  CO2  again 
proceeds  primarily  via  the  CO  oxidation  reaction  and  once  the  OH  radical  pool  increases 
after  the  disappearance  of  the  methane  the  CO  is  then  oxidized  to  form  CO2.  The  carbon 
dioxide  forms  downstream  of  the  regions  of  high  Ci/3,  and  CO  concentration. 

In  the  confined  flame  the  acceleration  of  the  gases  along  the  symmetric  axis  helped 
to  produce  a  flame  several  cm  in  length.  Without  the  additional  acceleration  near  the 
axis  of  symmetry  we  obtain  a  much  shorter  flame.  We  compute  a  flame  with  a  height 
of  approximately  1.25  cm.  Aside  from  the  height  of  the  two  flames,  the  most  dramatic 
difference  between  the  confined  and  unconfined  flames  is  in  the  fluid  dynamic  solution 
fields.  In  Figure  21  we  illustrate  the  stream  function  for  the  unconfined  flame.  We  observe 
that  the  streamlines  indicate  a  bowing  out  of  the  flame  above  the  inlet  with  a  gradual 
movement  of  the  flow  towards  the  symmetric  axis.  No  recirculation  cell  occurs  in  this 
configuration.  Also,  we  observe  in  Figure  22  that  there  are  two  large  vorticity  cells  that 
begin  approximately  4.0  cm  above  the  inlet  and  are  centered  approximately  0.5  cm  from 
the  axis  of  symmetry.  These  are  the  regions  of  highest  counterclockwise  rotation  in  the 
flame.  Here  the  flow  has  a  strong  velocity  component  towards  the  axis  of  symmetry.  This  is 
in  contrast  to  the  region  directly  above  the  inlet  where  the  vorticity  is  negative  and  the  flow 
direction  is  such  that  there  is  a  substantial  radial  component  of  the  velocity  causing  the 
flame  to  “bow  out” .  The  negative  vorticity  region  is  limited  to  the  area  directly  above  the 
jets.  At  a  height  of  approximately  1.0  mm  above  the  inlet,  the  vorticity  changes  sign  and 
the  inward  radial  component  of  the  velocity  begins  to  increase.  In  the  confined  flame,  due 
to  the  recirculation  cells,  the  high  vorticity  region  is  located  much  closer  to  the  symmetric 
axis. 

4.2.  Turbulent  Flames 

The  final  objective  of  our  Phase  II  work  plan  is  the  incorporation  of  a  k  -  c  turbulence 
model  in  our  fluid  dynamic-thermochemistry  equations.  The  turbulent  modeling  capabil¬ 
ities  of  AXIJET  lag  behind  the  laminar  modeling  capabilities  in  terms  of  physical  realism 
and  robustness  of  the  computational  algorithm,  even  though  they  are  close  to  the  state 
of  the  engineering  art.  This  is  because:  (a)  turbulence  is  inherently  a  three-dimensional 
time-dependent  phenomenon,  with  multiple  space  and  time  scales  that  are  likely  to  remain 
computationally  unresolvable  for  many  years,  even  with  optimistic  extrapolations  of  com¬ 
puting  power,  and  (b)  two-dimensional  time-averaged  models  of  turbulent  flows  contain 
extra  forms  of  coupling  not  found  in  laminar  flow  models  which  force  Newton’s  method 
to  proceed  at  such  small  implicit  time-steps  that  less  implicitly  coupled  techniques,  when 
tameable,  seem  preferable,  to  date,  in  terms  of  computing  resources. 
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A  simple  reminder  of  the  inaccessibility  of  turbulent  flow  computations  from  first 
principles  comes  from  a  spatial  scale  argument.  Let  L  be  an  integral  length  scale  for 
a  turbulent  flow,  say  the  diameter  of  a  pipe  bounding  the  flow,  and  let  l  be  the  scale 
of  the  smallest  turbulent  eddies  requiring  resolution  (the  Kolmogorov  scale)  in  a  direct 
Navier-Stokes  simulation  of  the  flow.  Then  L/l  =  0(Re3/4),  where  Re  is  a  Reynolds 
number  based  on  the  RMS  of  the  fluctuating  velocity  components  and  the  integral  length 
scale.  This  can  be  related  to  a  conventional  pipe  flow  Reynolds  number  upon  division 
by  an  average  turbulence  intensity,  typically  1  to  10  percent.  The  number  of  mesh  points 
required  is  0((L/f)3)  =  0(Re9/4).  For  a  Reynolds  number  of  104,  this  amounts  to  a  billion 
grid  points.  The  only  satisfactory  direct  simulations  of  turbulent  flows  performed  to  date 
have  been  for  Reynolds  numbers  of  O(l02). 

4.3.1.  Problem  Formulation 

Given  the  general  inadequacy  of  computational  methods  in  the  face  of  turbulence, 
it  is  preferable  to  modify  the  governing  equation  model  to  account  in  an  approximate 
manner  for  the  enhanced  mixing  that  occurs  in  turbulent  reacting  flows,  and  its  effects  on 
the  reaction,  than  to  neglect  turbulence  completely.  AXIJET-T  incorporates  a  standard 
Favre-averaged  k  —  c  flow  model  and  a  standard  mixture-fraction-covariance  model  for  the 
chemical  reaction.  The  k  —  «  model  reduces  the  evaluation  of  the  turbulent  stresses  to  a 
pair  of  fields  which  satisfy  transport  equations  similar  to  those  governing  steady  laminar 
flows.  These  are  combined  to  generate  local  time  and  length  scales  for  the  turbulence. 
The  mixture-fraction-covariance  model  reduces  the  evaluation  of  the  local  thermodynamic 
state  of  the  flow  to  another  pair  of  fields,  denoted  /  and  g  herein,  which  also  satisfy  steady 
transport  equations.  Within  AXIJET-T’s  k-t~f-g  formulation  there  are  three  main  sub¬ 
hypotheses:  an  isotropic  eddy  viscosity,  a  /9-function  probability  density  function  for  the 
mixing  of  the  fuel  and  oxidizer  streams,  and  the  laminar  flamelet  hypothesis.  Discussion 
of  these  submodels  is  postponed  until  the  following  section. 

Since  AXIJET-T  represents  the  flow  field  in  terms  of  steady,  time-averaged  quanti¬ 
ties,  precise  definition  of  the  averaging  employed  is  in  order.  For  constant  density  flows, 
turbulence  modeling  is  based  on  the  decomposition  introduced  by  O.  Reynolds  (1895)  of 
a  time-dependent  field  <f>(x,t)  into  a  steady  time-averaged  part  (denoted  with  an  overbar) 
and  a  fluctuating  part  with  zero  mean  (denoted  with  a  prime): 

<^{x,f)  =  <?(x)  +  <£'(x,f) 

where 

T 

#(*)  =  lim  i  f  i{x,t)  dt 
r-oo  i  Jo 

For  variable  density  flows,  Reynolds  averaging  of  the  governing  equations  is  inconve¬ 
nient  in  that  many  correlations  between  density  and  other  fields  emerge.  In  combustion 
applications,  the  density  may  vary  by  an  order  of  magnitude  over  the  flow  domain,  so  that 
these  terms  are  non-negligible.  The  density-weighted  time  averaging  introduced  by  Favre 
(1969)  (53)  overcomes  this  inconvenience  for  variable  density  flows  (although  it  still  requires 
modeling  of  triple  correlations  in  the  turbulent  stress  terms).  The  Favre  decomposition 


is  written  in  terms  of  a  density-weighted  average  (denoted  with  a  tilde)  and  a  fluctuating 
part  with  zero  density- weighted  mean  (denoted  with  a  double  prime): 

^(x,t)  =  4>{x)  +  t) 

where 

T 

^(x)  =  ^  X  Jim  i  f  P<t>{x,  t)  dt 
P  T—oc  1  J  o 

The  following  simple  relationship  between  Reynolds  and  Favre  averages  shows  that 
they  reduce  to  the  same  decomposition  when  density  is  constant: 

4>  =  $(x)  +  p><VIp 

In  terms  of  Favre  averages,  the  primitive  variable  governing  equations  in  two-dimensional 
axisymmetry,  and  in  the  subsonic  regime  for  which  AXIJET  is  designed,  are  given  by: 
Continuity: 

a  .  .  i  a  , 

-(/.«)  +  7rrM~o 

Axial  Momentum: 

a  __  — ——  i  a .  __  — dp  dfxx  i  a ,  . 

—  (/Juu  +  pu"u ")  +  -^P™  +  rpv"u")  +  ~  +  -^(rf,r)  =  0 

Redial  Momentum: 


a  .  __  1  3.  — jj-rj.  dp  Id.  .  Tgg  dfra 

Tz(PUV  +  PU  V  )  +  rrr{,PVV  +  ,PV  V  )  +  Tr+  7fr[rTrr)  “  T  +  ~57  =  ° 

No  gravity  body  force  term  is  retained  in  the  momentum  equations  because  the  large 
forced  convective  flow  velocities  and  short  convective  time  scales  typical  of  turbulent  ap¬ 
plications  render  it  insignificant. 

The  time-averaged  molecular  stress  terms  are  modeled  by  the  usual  Newtonian 
formulae  in  terms  of  time-averaged  velocity  gradients: 


,  aa  2 

r**  =  -M-r*  ~  °i 

{dV  2  . 

rrr  =  -M(2_-5V.,1| 

T)t  =  -  |v  •«) 


.da  dv, 

r”  = -  -"I  Tr  +  5I1 
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The  triple-correlation  turbulent  stress  terms  are  modeled  by  the  usual  Doussinesq 
hypothesis  in  analogy  to  the  molecular  stresses  above: 


/>«"■*"  =  +  \tk 

+  5  Pk 

- rj—r.  ,du  dv. 

=  -»>\Tr  +  Tz  1 

The  turbulent  stress  terms  introduce  a  turbulent  (or  eddy)  viscosity  fit  in  analogy  with 
the  laminar  (or  molecular)  viscosity  fi.  We  note  that  the  a  and  0  required  in  the  molecular 
stress  terms  are  not  available  without  additional  modeling  assumptions.  However,  the 
molecular  stresses  will  be  generally  negligible  in  comparison  with  the  turbulent  stresses 
except  near  walls,  since  ft  «  fit.  Therefore,  we  may  use  u  and  v  in  place  of  a  and  0 
in  the  molecular  terms  and  sum  the  two  stresses  away  from  walls.  In  wall  regions,  the 
small  errors  thus  committed  are  acceptable  in  view  of  our  reliance  upon  standard  constant 
density  wall  function  approximations,  themselves. 

In  terms  of  the  Favre  velocity  fluctuations,  we  may  define  a  specific  turbulent  kinetic 
energy 

k  =  ^(ti"2  +  v"2)  (units:  l2t~2 j 

m 

and  a  specific  turbulent  kinetic  energy  dissipation  rate 


e  =  v  < 


(du"\2  s 

(  dv"\2  (% 

r"\2 

(  dv 11 

du"\2\ 

2 

[U-j + 1 

U7)  +(" 

r) 

+  (a7 

+  ar)  ) 

[units:  l2t~3\ 


The  eddy  viscosity  will  be  specified  in  terms  of  k  and  e,  so  the  non-reacting  turbulence 
model  is  closed  by  the  following  two  transport  equations: 

Turbulent  Kinetic  Energy: 


d  ,  ...  1  3  .  ...  d 


l±f  f*e  dk\ 
rdr  V  okdr) 


=  $  —  Pf. 


Turbulent  Kinetic  Energy  Dissipation  Rate: 


d  .  _  ,  ia. 

+  rSirpve) 


r  dr  \  otdr  ) 


These  equations  require  definition  of  diffusive  transport  coefficients  in  terms  of  the 
effective  viscosity 

fie  =  fi  +  fit 

and  k  and  e  “Prandtl”  numbers  ok  and  o>,  which  we  take  as  standard  fits  from  constant 
density  wall-bounded  shear  fiow  experience  with  the  model: 
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<7*  =  1.0 
a(  =  1.3 

Each  equation  has  a  source  term  based  on  the  local  dissipation  of  the  mean  flow: 


— {•[(£)'*  (£)’* (')']*  (S ■ * £)’}  i— 

and  a  corresponding  sink  term.  Whereas  the  derivation  of  the  ^-equation  rests  on  good 
theoretical  grounds  (at  least  for  isotropic  turbulence),  the  derivation  of  the  {-equation  is 
weaker  and  relies  on  dimensional  analysis  and  empirical  fitting.  The  constants  we  employ 
in  the  source  and  sink  term  are  again  standard  fits  from  the  literature: 

Cel  =  1.47 
Ct2  =  1.92 

Note  that  ratio  k/e  provides  a  characteristic  temporal  scale  for  the  turbulence  and  fc3/2/e 
a  characteristic  spatial  scale. 

We  model  the  chemical  composition  of  the  reacting  flow  by  means  of  a  Favre-averaged 
mixture  fraction,  /,  where  /(x)  is  defined  as  the  fraction  of  the  mixture  at  x  originating 
from  the  fuel  stream,  and  hence  varies  from  a  value  of  0  at  the  oxidizer  inlet  to  1  at  the 
fuel  inlet.  An  unambiguous  definition  of  /  can  be  made  in  terms  of  the  local  molecular 
composition,  the  molecular  weight,  and  the  atomic  composition  of  each  molecule  involved 
in  the  reaction  mechanism.  For  concreteness,  an  example  based  on  a  pure  methane  fuel 
stream  and  an  oxygen/nitrogen  oxidizer  stream  follows.  Let  there  be  K  molecular  species 
in  a  reaction  mechanism,  with  molecular  masses  M*  and  local  mass  fractions  Y*.  Let  a *  q 

denote  the  number  of  carbon  atoms  in  the  It4*1  molecule  and  a*  jy  the  number  of  hydrogen 
atoms.  Finally,  let  Mq  and  Mjj  denote  the  atomic  masses  of  carbon  and  hydrogen.  Then 

/(jt)  =  £.LA“k,cMcYk(x)/Mk)  + 

ych< 

Two  particular  values  of  /  axe  of  special  interest,  the  asymptotic  exit  plane  value 


f ezit  — 


1  +  mo/m  p 


and  the  stoichiometric  value 


fSt  - - - - 

1  +  {uqMq/upMp) 

where  mp  and  m0  are  the  integrated  mass  fluxes  of  fuel  and  oxidizer,  i/p  and  uq  their  sto¬ 
ichiometric  coefficients  (adjusted,  if  necessary,  for  the  fact  that  either  stream  may  contain 
a  nonzero  amount  of  inert),  and  Mp  and  Mq  their  molecular  masses. 
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Under  the  assumption  that  all  state  variables  of  the  flow  are  functions  of  /,  solution 
of  a  source-free  transport  equation  for  /  describes  the  complete  thermodynamic  state  of 
the  mixture.  However,  because  of  the  underlying  turbulence,  it  is  unreasonable  to  assume 
that  f(x)  is  a  deterministic  function,  and  it  is  typically  modeled  as  a  random  function 
governed  by  some  Favre  probability  density  function  P{f).  Therefore,  any  state  variable 
or  function  thereof  (such  as  temperature  or  viscosity)  is  derivable  as  a  probability  integral 
over  the  domain  of  allowable  values  of  /: 

Jo 

The  density  is  the  lone  exception  to  this  formula,  because  of  the  density-weighted 
averaging  employed.  The  appropriate  formula  is  instead: 

rL  l7U)p{fW 

Evaluation  of  P(f)  typically  requires  a  set  of  moments  of  /.  The  simplest  schemes 
involve  only  one  higher  moment,  the  Favre  covariance,  denoted  g  and  defined  by: 

g  =  (n 2 

where  f"  =  /  -  /.  The  specific  choice  of  P  enters  only  into  a  free-standing  module  bundled 
with  AXIJET-T,  which  in  turn  affects  the  main  AXIJET-T  code  only  through  a  data  file 
giving  the  state  functions  necessary  to  advance  the  calculation  in  terms  of  /  and  g. 

The  Favre-averaged  mixture  fraction  satisfies  the  following  homogeneous  equation: 
Mixture  Fraction: 


The  mixture  fraction  covariance  satisfies  a  modeled  equation  similar  in  structure  to 
the  e  equation: 

Covariance  of  Mixture  Fraction  (“Unmixedness”): 


d_ 

dz 


.  „  ,  la,  a 

(<>“»)  +  ;*<’•*'»)  -  Tt 


1  d_  f  /*e 
r  dr  \  og  dr  ) 


=  C, it,  -  C,2^- 


Ab  with  the  k  —  e  pair,  modeling  constants  are  required,  two  “Prandtl  numbers”: 

Of  =  0.6 
Og  —  0.6 
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and  source  and  sink  term  coefficients,  for  which  we  employ  the  fits: 


Cg  1  =  2.8 

Cg2  =  2.0 

The  source  term  consists  of  squares  of  gradients  of  /,  in  analogy  to  the  mean  dissipation 
term  of  the  k  —  e  model: 


$g  = 


[units:  ml  3t  *] 


4.2.2.  Transport  and  Chemistry  Model 

The  transport  and  chemistry  model  consists  of  three  independent  models:  a  set  of 
deterministic  state  functions  of  mixture  fraction,  a  probability  density  function  with  which 
these  are  weighted  in  constructing  the  mean  properties,  and  a  formula  for  the  enhanced 
mixing  due  to  the  turbulence  in  the  form  of  nt‘ 

For  the  first,  an  assumption  of  thermo-chemical  equilbrium  is  often  asumed  in  the 
literature,  as  matter  of  convenience  in  view  of  the  ignorance  of  finite-rate  kinetics.  Though 
an  equilibrium-based  data  file  can  be  prepared  for  the  methane-air  flame  demonstrated  in 
the  results  section  below,  we  employ  instead  a  primitive  laminar  flamelet  hypothesis,  which 
takes  advantage  of  our  earlier  work  in  the  modeling  of  finite-rate  methane-air  flames.  The 
laminar  flamelet  model  is  based  on  a  picture  of  a  turbulent  flame  front  as  a  convoluted 
ensemble  of  laminar  flamelets  -  counterflowing  diffusion  flames,  to  be  precise.  The  local 
composition  is  taken  to  the  composition  of  the  corresponding  diffusion  flame  at  the  same 
value  of  the  mixture  fraction.  Refinements  to  the  laminar  flamelet  model  take  into  account 
variations  in  the  local  fluid  rate  of  strain.  This  capability  is  not  incorporated  into  the 
present  version  of  AXIJET-T. 

Because  of  the  way  AXIJET-T  is  formulated  above,  knowledge  of  most  state  variables 
is  not  needed  to  advance  the  calculation;  only  density  and  viscosity  are  needed  at  each 
intermediate  iteration.  Consequently,  the  evaluation  of  the  temperature  and  mass  fractions 
of  various  species  can  be  deferred  to  a  post-processing  step. 

Plots  of  density  and  viscosity  as  a  function  of  mixture  fraction  are  given  in  Figures  23 
and  24  following  and  a  plot  showing  the  temperature  and  the  two  major  reactants  is  given 
in  Figure  25.  The  data  comes  from  the  Tsuji  flame  configuration  solved  in  [54].  Note  that 
for  this  particular  flame,  the  stoichiometric  value  of  /  is  0.0551,  and  this  is  where  the  peak 
of  the  temperature  profile  (also  the  trough  of  the  density  profile)  sits. 

For  the  second  model,  namely  the  probability  density  function  for  /,  “clipped  Gaus- 
sians"  or  Incomplete  Beta  functions  are  often  used.  The  latter  is  built  into  the  AXIJET 
preprocessor.  Its  definition  is 


-/)•-*  <J/’ 


(4.92) 
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where  the  exponential  parameters  are  derived  from  the  Favre  mean  and  covariance: 


a  =  / 


/(l  -  /) 


-  1 


.  *=(!-/) 


/(I  -  /) 


9 


-  1 


(4.93) 


It  is  readily  verified  for  o  >  0  and  6  >  0  that  P  satisfies  the  basic  properties  of 
a  probability  density  function  over  the  interval  (0, 1):  it  is  non-negative  and  has  a  unit 
integral  over  the  complete  interval.  Note  that  g  must  be  confined  to  the  range  0  <  g  < 
/( 1  —  /)  for  this  to  be  true.  Very  small  values  of  g  lead  to  delta-function-like  peaked 
distributions  at  the  endpoints,  and  g  =  0  represents  certainty  that  the  flow  is  in  the  state 
characterized  by  /.  Large  values  of  g  lead  to  Gaussian-like  distributions  about  /.  Some 
sample  plots  of  P  for  different  /  and  g  are  contained  in  Appendix  A.  Appendix  A  also 
addresses  some  fine  points  in  the  numerical  quadrature  of  integrals  containing  /^-functions. 

We  follow  standard  practice  in  adopting  for  our  third  model  an  isotropic  turbulent 
viscosity  coefficient 

Turbulent  Viscosity: 


Ht  =  CpPk2/*  [units:  ml  J] 
which  contains  another  fitted  parameter: 

Cp  =  0.09 

4.3.3.  Boundary  Conditions  and  Wall  Functions 

Discussion  in  this  section  is  patterned  after  the  geometry  that  is  presently  hard-coded 
into  the  release  copy  of  AXIJET-T,  namely  that  of  a  two-fluid  jet  with  a  coaxial  coplanar 
inlet  boundary  upstream,  adiabatic  no-slip  bounding  walls,  and  an  outflow  boundary. 
(Figure  26  applies  equally  well  in  this  section.)  An  ( z,r )  coordinate  system  is  adopted 
with  its  origin  at  the  center  of  the  inlet  plane  and  key  dimensions  as  follows.  The  radius 
of  the  cylindrical  fuel  jet  is  rjn  and  the  outer  radius  (which  is  the  same  as  that  of  the 
annular  oxidizer  jet)  is  rout  The  outflow  is  at  zout.  In  describing  phenomena  in  the  wall 
region,  a  local  coordinate  system  is  also  adopted  at  wall  boundaries  in  which  y  is  a  normal 
coordinate  with  its  origin  at  the  wall  itself,  and  and  U||  are  the  local  normal  and 
tangential  velocities,  respectively. 

Due  to  radial  symmetry,  only  half  of  an  (z,r)  plane  need  be  considered.  Thus  four 
types  of  boundaries  are  discussed  in  conjunction  with  the  particular  geometry  currently 
in  AXIJET-T:  inflow,  outflow,  symmetry  and  wall.  These  are  the  same  types  likely  to 
be  required  in  most  future  applications  of  AXIJET-T,  except  possibly  for  permeable  wall 
boundaries  (along  which  uy  =  0,  but  «_l  ^  0)  and  nonadiabatic  wall  boundaries.  We  begin 
with  wall  boundaries,  because  they  are  both  key  and  problematic  features  in  turbulence 
modeling,  and  because  results  from  the  wall  boundaries  are  subsequently  used  at  inlet 
boundaries  a a  well.  From  this  subsection  onward,  ail  tilde  notation  of  the  Favre  averaging 
is  dropped  but  implicitly  understood. 
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4. 2.3.1.  Wall  Boundaries 

Wall  boundary  conditions  along  an  impermeable  no-siip  surface  are  straightforward, 
namely  U||  =  =  0,  whence  k  =  t  =  0.  If  the  wall  is  adiabatic,  we  also  have  that 

yL  —  o  and  g  =  0.  In  the  history  of  the  numerical  simulation  of  high  Reynolds  number 
turbulent  flows  it  has  often  difficult  for  economic  reasons  to  resolve  the  wall  boundary 
layer  sufficiently  well  for  these  straightforward  boundary  conditions  to  be  applied  without 
severely  contaminating  the  accuracy  of  the  solution  outside  of  the  layer.  Furthermore,  the 
process  of  resolving  the  wall  region  accurately  enough  to  capture  the  normal  gradients, 
inevitably  introduces  large  numbers  of  points  in  the  low  Reynolds  number  region  near  the 
wall,  which  increases  the  volume  of  computational  work  required.  Furthermore,  in  the  low 
Reynolds  number  regions,  the  k  —  c  model,  and  the  resulting  estimate  for  fit  are  invalid, 
and  must  be  replaced  with  additional  local  assumptions  anyway. 

Therefore,  it  is  customary  to  look  to  empirical  fits  of  experimental  data  to  interpolate 
between  the  bulk  high  Reynolds  number  turbulent  flow  and  the  no-slip  wall.  This  practice 
is  complicated  by  the  fact  that  the  bulk  flow  depends  on  the  field  variables  at  the  outer  edge 
of  the  wall  function  region,  and  the  parameters  of  the  wall  functions  are  defined  in  terms 
of  the  bulk  flow  -  the  patching  process  must  be  implicitly  handled.  (This  requirement 
for  implicitness  is  no  different  a  situation  than  that  encountered  in  enforcing  a  system  of 
elliptic  PDEs  without  wall  functions  all  the  way  to  the  wall,  but  it  complicates  the  coding 
in  that  a  transition  region  in  which  both  sets  of  formulae  apply  must  be  determined  in 
order  for  the  solution  process  to  go  forward.) 

The  most  convenient  links  by  which  the  patching  process  is  mediated  are  the  shear 
stress  at  the  wall,  denoted  herein  by  rWt  and  the  “friction  velocity" ,  denoted  Ur  and  defined 
in  terms  of  rw  by  Ur  =  \ZrU)/pw.  These  are  used  to  define  a  host  of  nondimensional 
quantities  as  follows: 


The  first  of  these  quantities,  y+,  is  used  to  divide  the  domain  into  three  layers,  ac¬ 
cording  to  classical  turbulent  flow  taxonomy  (the  transitions  between  these  layers  are  ap¬ 
proximate.)  The  laminar  (or  Couette)  layer  is  0  <  y+  <  12;  the  log  layer  is  12  <  y+  <  300; 
and  the  outer  layer  is  y+  >  300.  The  values  of  k  and  e  at  grid  points  at  and  near  the  wall 
will  depend  on  the  location  of  the  grid  points  in  terms  of  y+.  (We  note  that  the  physical 
locations  (in  y)  a  fixed  y+  change  as  the  flow  evolves,  since  y+  depends  on  the  local  shear 
stress  at  the  wall,  which  varies  with  bulk  flow  conditions.)  The  “law” 


in  which  k  =  0.42  and  B  =  5.5,  holds  in  the  log  layer  and  constitutes  a  transcendental 
equation  for  UT ,  which  may  be  solved  by  iteration  at  any  point  in  the  layer.  Having  deter¬ 
mined  the  dimensionless  normal  distance,  we  are  able  to  complete  the  following  scheme: 
Parallel  Velocity: 


u+  =  /  v+. 

\  In  v+/k  -f  B , 


Specific  Turbulent  Kinetic  Energy: 


if  y+  <  12  1 
if  y+  >  12  / 


0.05(y+)2, 

1.25  +  0.325(y+  -  5), 
4.5  -  (y+  -  15)/37.5, 
3.3, 


if  y+  <  5  \ 

if  5  <  y+  <  15  I 
if  15  <  y+  <  60  [ 

if  60  <  y+  <  200  J 


Specific  Dissipation  Rate  of  Turbulent  Kinetic  Energy: 


(4.94) 


(4.95) 


+  [  0.1  4-  y+/120,  if  y+  <  12  1 

\  l/(«y+),  if  y+  >  12  / 


(4.96) 


The  log  law  is  illustrated  in  Figure  27  and  these  fits  for  fc+  and  c+  in  Figures  28  and 
29. 

One  final  feature  of  the  wall  function  technique  is  the  Van  Driest  formula  for  fit  within 
the  log  layer: 

Wall  Region  Turbulent  Viscosity: 


fit  =  £(\/l  +  4a2  -  1)  (4.97) 

it 

where 

a  =  *y+(  1  -  e~v+/A+)  , 

in  which  A+  =  26.7. 

4. 2.3.2.  Inflow  Boundaries 

Generally,  normal  and  tangential  mean  velocities,  as  well  as  density  are  available  at 
inflow  boundaries.  In  multicomponent  problems,  the  mixture  fraction  /  is  often  known 
with  precision  in  inlet  streams,  in  which  case  its  covariance,  g,  is  zero.  In  the  absence  of 
measured  data,  specific  turbulence  kinetic  energy  can  be  be  hypothesized  to  be  a  fraction  of 
the  specific  mean  kinetic  energy,  k  =  c(ujj  +  u2  )  (c  is  dimensionless).  Again,  in  the  absence 

of  superior  information,  a  constant  turbulence  length  scale  hypothesis  yields  c  =  it3/2//  (/ 
has  dimensions  of  length).  We  use  these  formulae  towards  the  middle  of  the  inlets,  as 
described  below,  assuming  that  the  inlet  profiles  are  themselves  well  developed.  We  choose 
the  constants  by  grafting  them  onto  boundary  layer  formula  at  the  inlet  edges.  Hence,  all 
inflow  data  are  Dirichlet. 

In  the  case  at  hand,  the  inlets  are  z  =  0,  0  <  r  <  rjn  (cylindrical  fuel  inlet,  subscript 
F)  and  z  —  0,  rjn  <  r  <  rout  (annular  oxidizer  inlet,  subscript  O). 

For  fully-developed  turbulent  flow,  we  have,  away  from  the  wall  regions: 
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Velocities: 


vr  —  0 


v,  =  VF,mag(--V-)1/7  ,  0  <  r  <  rj 


v*  —  vO,max( 


T“ - )^?  »  rin  <  r  <  ^rout 

''f’out  —  rin 


vz  —  vO,max( - T - *  ^rout  <  r  <  rout 

r0ut  —  ''f’out 

where  the  maxima  at  r  =  0  and  r  =  Ar0„t  are  given  by 

,60, 

VF,max  =  (49^°^ 

,8°,  1 

v°,max  -  (49'  1  +  A/(7(l  +  rin/r0„t))  0 

where  Op  and  Vq  are  the  fuel  and  oxizider  mass  fluxes  divided  by  the  inlet  densities  and 
areas,  and  A  =  J 


areas,  and  A  = 

By  equating  these  expressions  for  the  parallel  velocity  components  away  from  the  inlet 
walls  with  log  law  velocity  profiles  at  a  given  matching  value  of  y+,  we  may  deduce  the  local 
friction  velocity  Ur  at  each  of  the  three  inlet  walls,  r  slightly  less  than  in  the  fuel  jet, 
r  slightly  greater  than  rin  in  the  oxidizer  jet,  and  r  slightly  less  than  rout  in  the  oxidizer 
jet.  From  these  values  of  UT,  the  axial  velocity  profiles  can  be  corrected  near  the  wall 
through  (4.94).  A  value  of  about  200  for  y+  is  reasonable  for  this  matching.  Because  of 
the  wall  correction  to  the  l/7th  power  law  velocity  distribution,  the  mass  flux  through  each 
inlet  will  be  slightly  reduced  from  the  design  mass  fluxes,  from  which  wo.ma*  an<*  VF,max 
are  initially  calculated  under  the  assumption  of  a  full  profile.  It  is  therefore  advisable  to 
renormalize  the  velocity  distribution  once  it  is  in  final  shape.  (In  typical  applications,  this 
effect  is  2%  or  less,  and  the  process  is  not  iterated.) 

Since  UT  is  known,  k ,  e  and  m  can  be  set  in  each  wall  inlet  layer  by  means  of  (4.95) 
through  (4.97).  This  process  provides  three  sets  of  values  of  k  and  e  at  y+  =  200  which 
can  then  be  used  to  determine  the  constants  c  and  l  in  the  sets  of  formulae 
Specific  Turbulent  Kinetic  Energy: 

k  =  cVg 

Specific  Rate  of  Dissipation  of  Turbulent  Kinetic  Energy: 


Of  course,  since  the  annular  inlet  provides  a  matching  at  both  sides,  the  constants  c 
and  l  may  be  slightly  different  when  fit  independently  from  either  side.  In  practice,  this 
discrepency  should  be  modest  but  it  is,  in  any  event,  dispensed  with  by  calculating  the  k 
and  c  distributions  from  both  ends  and  blending  them  linearly  in  the  middle. 
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Finally  we  have  the  inlet  conditions  for  the  f  —  g  submodel: 
Mixture  Fraction: 


1, 

1  fsti 

0. 


if  0  <  r  <  rjn 
if  r  =  rin 
if  rjn  <  r  <  rou  t  , 


Covariance  of  Mixture  Fraction: 


<7  =  0 


(4.98) 


4. 2.3.3.  Outflow  Boundaries 

The  outflow  boundary  is  assumed  to  be  far  enough  downstream  that  zero  gradient 
conditions  can  be  imposed  on  the  radial  momentum,  and  all  other  field  variables  except 
the  axial  momentum.  The  flow  is  overconstrained  by  a  zero  gradient  condition  on  the  axial 
momentum  which  translates  through  the  continuity  equation  and  the  lateral  boundary 
conditions  to  a  statement  that  vr  =  0.  Therefore,  the  axial  momentum  is  required  to 
satisfy  the  less  stringent  condition  of  zero  curvature.  In  the  case  at  hand,  the  outlet  is 
z  =  *ont»  0  <  r  <  rout.  We  have: 

Velocities: 


d(/™r) 

dz 


All  other  variables  <f>: 


d2(pv«) 

dz* 


=  0 


dj 

dz 


=  0 


4. 3.3.4.  Symmetry  Boundaries 

On  a  symmetry  boundary  the  fields  satisfy  zero-normal  gradient  conditions.  In  ad¬ 
dition,  the  velocity  normal  to  the  symmetry  boundary  must  be  zero  by  continuity.  We 
have: 

Normal  velocity,  vr: 

vr  =  0 


All  other  variables  <f> : 


d4> 

dr 


=  0 
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4.2.4.  Initial  Conditions 

It  is  useful  to  calculate  outflow  parameters  based  on  the  assumption  of  fully  developed 
turbulent  outflow  profiles,  since  this  provides  information  which  can  be  used  in  constructing 
an  initial  iterate  by  axial  interpolation  of  inlet  and  outlet  streams.  For  this  purpose,  we 
assume  a  fully  developed  outflow  and  construct  an  axial  velocity  in  direct  geometric  analogy 
to  the  fuel  inlet,  except  that  the  radius  is  different  and  the  mass  flux  is,  of  course,  the  sum 
of  that  of  both  inlet  streams.  As  at  the  inlet,  we  assume  that  g  =  0,  and  we  set  /  =  /eant 
at  all  r.  All  fields  except  for  k  and  c  are  then  linearly  interpolated  in  z  between  inlet  and 
exit.  Consulting  the  state  relationship  of  p(f)  provides  a  global  estimate  for  the  density. 
It  is  then  possible  to  construct  an  estimate  for  the  axial  momentum  pvt  at  all  points,  and 
thereby  of  the  integrated  axial  mass  flux  at  each  station.  Since  these  mass  fluxes  by  station 
must  match  the  sum  of  the  inlet  fluxes,  the  velocity  profiles  are  normalized  by  the  ratio  of 
fluxes  in  an  approximate  attempt  to  begin  the  calculation  with  a  flow  field  which  is  axially 
non-divergent. 

The  turbulence  parameters  are  set  at  their  lower  bounds  everywhere  except  at  the 
inlets  in  the  initial  estimate,  so  that  the  flow  effectively  begins  as  a  laminar  one.  Lin¬ 
ear  interpolation  between  assumed  turbulence  profiles  at  the  inlet  and  outlet  stations 
was  attempted  and  found  highly  unsatisfactory  because  in  typical  applications  turbulent 
"hot  spots”  tend  to  be  localized  in  the  axial  coordinate,  are  thus  poorly  approximated 
by  a  monotonic  interpolation  between  the  inlets.  Furthermore,  their  appearance  in  sev¬ 
eral  “guessed”  source/sink  terms  pulls  the  flow  in  several  directions  at  the  outset,  which 
confuses  the  solution  algorithm  in  its  initial  pressure  equilibration  progress. 

4.2.5.  Computational  Approach 

AXIJET-T  is  basically  a  SIMPLER  (Semi-Implicit  Method  for  Pressure-Linked  Equa¬ 
tions  -  Revised)  algorithm  after  Patankar  (55]  based  on  a  staggered-grid  Finite  Analytic 
discretization  of  the  governing  equations  after  Chen  and  Chen  [56]  (see  also  [57]  for  a 
turbulent  flow  application).  The  major  refinement  over  a  conventional  SIMPLER  imple¬ 
mentation  is  the  use  of  a  direct  solver,  YSMP  (Yale  Sparse  Matrix  Package  [58])  to  handle 
the  pressure  and  pressure  correction  equations  without  iteration.  We  repeat  that  all  tilde 
notation  of  the  Favre  averaging  is  dropped,  since  the  algorithm  applies  equally  well  to  any 
physical  formulation  which  fits  into  the  generic  five-point  transport  operator  on  a  staggered 
grid. 

As  applied  to  the  extended  system  of  field  equations  (for  u(=  va ),  t/(=  t/r),  p,  k ,  c, 
/,  and  g ),  and  algebraic  relationships  (for  p,  p,  and  /*;),  the  SIMPLER  algorithm  has 
the  following  form  (consult  [55]  for  the  well-documented  details).  The  roles  of  the  lin¬ 
earized  velocity  components  u  and  v  at  each  iteration  are  played  by  two  sets  of  fields: 
“pseudo-velocities”  (6  and  6)  which  satisfy  the  discrete  momentum  equations  without  the 
pressure  gradient  terms,  and  “true”  velocities  which  satisfy  the  full  momentum  equations, 
albeit  based  on  a  provisional  pressure.  Similarly,  the  pressure  is  decomposed  into  two 
fields:  a  provisional  pressure  (p*)  which  satisfies  the  Poisson  equation  obtained  by  taking 
the  divergence  of  the  momentum  equations  with  the  source  term  computed  based  on  the 
pseudo-velocities  only,  and  a  pressure  correction  ( p ')  based  on  the  same  discrete  Poisson  op¬ 
erator,  but  with  the  “true”  velocities  used  in  computing  the  source  term.  As  implemented 
in  AXIJET-T,  it  takes  the  following  form: 
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The  “SIMPLER”  Algorithm 


(0)  Specify  ufc,  vk,  (p*)*,  etc.,  at  iteration  index  k  =  0 

(1)  Solve  linearized  momentum  eqns.  for  temporary  velocities  u*,  v* 

(2)  Solve  pressure-correction  eqn.  for  (p/)*+1 

(3)  Solve  velocity-correction  eqns.  for  using  (p,)*+1 

(4)  Solve  pressure-free  momentum  eqns.  for  <1,  0 

(5)  Solve  for  provisional  pressure  (p*)*+1 

(6)  Solve  linearized  conservation  eqns.  for  other  fields 

(7)  Update  state  and  transport  properties 

(8)  If  converged,  STOP 

(9)  k  *-  k  +  1;  go  to  (l) 


The  cycle  (1-9)  is  ordered  differently  from  the  original  version  of  SIMPLER  in  two 
ways.  In  the  original,  Steps  (6)  and  (7)  immediately  follow  step  (3).  The  two  orderings  are 
nearly  equivalent,  however,  since  the  only  fields  modified  between  steps  (3)  and  (6)  are  the 
pseudo-velocities  and  the  pressure,  on  which  none  of  the  other  field  variables  depend,  and 
the  dependence  of  the  coefficients  of  the  equations  used  in  steps  (4)  and  (5)  upon  the  other 
fields  may  be  suppressed  anyway,  to  avoid  recomputation.  Apart  from  the  unimportant 
difference  in  the  cyclic  order,  the  cycle  begins  at  a  different  stage  in  the  original,  namely 
at  our  step  (4). 

Each  of  stages  (1)  through  (6)  involves  solving  large  sparse  linear  systems  of  algebraic 
equations  for  the  discrete  unknowns.  The  equations  arising  in  (2)  and  (4)  are  solved  by 
direct  Gauss  elimination,  which  incurs  a  large  storage  overhead,  but  should  be  regarded  as 
a  minimum  step  in  the  direction  of  implicitness  for  SIMPLE-type  algorithms  on  large  grids, 
over  which  relaxation  methods  are  too  inefficient  to  be  worth  the  storage  savings.  The  other 
systems  are  solved  by  block  line  Gauss-Seidel,  sweeping  in  the  convective  direction.  Sweeps 
are  repeated  until  the  residuals  of  the  linearized  discrete  equations  for  each  field  in  turn 
satisfy  an  absolute  convergence  criterion.  This  criterion  is  set  stringently  in  the  release 
version  of  AXIJET,  since  typically  only  a  handful  of  sweeps  are  necessary  in  combination 
with  the  below-mentioned  implicit  time-stepping.  Introduced  for  nonlinear  stability,  it 
has  the  side-effect  of  providing  diagonal  dominance  which  enhances  the  efficiency  of  the 
relaxation  sweeps. 

In  an  effort  to  control  the  notorious  instability  of  the  SIMPLER  class  of  algorithms  in 
the  presence  of  a  variable  density  k  —  c  model,  AXIJET  incorporates  sink-term-iinearization 
in  the  k,  e  and  g  equations,  and  implicit  time-stepping  in  ail  equations  except  for  the 
Poisson  pressure  equation.  In  addition,  “filtering”  is  performed,  if  necessary,  once  each 
cycle  through  the  SIMPLER  loop  to:  (a)  bound  the  k  and  e  fields  away  from  zero  or 
negative  values,  by  clipping  them  into  the  interval  [lO-10,oo)  (b)  clip  the  /  and  g  fields 
into  their  doubly  bounded  intervals  of  (0,  l]  and  [0,  /(I  —  /) j,  respectively.  While  doing 


the  filtering,  the  code  checks  that  the  outflow  boundary  is  truly  an  outflow  (i.e.,  that 
v*(*out)  >  0).  (A  warning  message  is  printed  if  this  condition  fails,  but  no  corrective 
action  is  taken,  since  it  is  likely  a  symptom  of  invalid  placement  of  an  outflow  boundary  if 
it  persists.) 

Dirichlet  boundary  conditions  arc  explicitly  incorporated  into  the  governing  equation 
set,  while  the  Neumann  symmetry  conditions  are  handled  by  solving  a  discrete  form  of  the 
governing  equations  on  the  axis  itself,  by  means  of  defining  a  row  of  auxiliary  boundary 
cells  “below”  the  symmetry  axis.  Neumann  outflow  conditions  are  incorporated  by  solving 
discrete  equations  up  to  the  penultimate  station  and  extrapolating  the  unknowns  of  the 
final  station. 

The  wall  functions  are  enforced  in  a  region  near  r  =  rout,  0  <  z  <  z0«t  by  explicitly 
overwriting  the  k ,  e,  nt  and  va  fields  in  the  last  few  discrete  rows.  (The  exact  number 
of  rows  is  a  user-specifiable  parameter).  The  transcendental  equation  for  UT  is  solved 
by  a  Newton-Rap hson  iteration  at  each  station  at  each  iteration.  The  trrnsverse  extent 
of  the  “swept”  region  solved  by  the  block-line  Gauss-Seidel  algorithm  is  trui  *  a  ed  away 
from  the  wall  region  for  the  k  and  e  fields,  since  the  governing  equations  from  which  the 
discretizations  are  derived  are  not  valid  in  this  low  Reynolds  number  region. 

Discretization  of  the  governing  partial  differential  equations  is  performed  with  the 
Finite  Analytic  method,  which  may  be  thought  of,  as  far  as  the  structure  of  its  dis¬ 
crete  operator  is  concerned,  as  an  interpolated-upwinding  Finite-Difference  scheme.  Like  a 
lowest-order  primitive  variable  Finite-Difference  scheme  for  this  set  of  governing  equations, 
it  makes  use  of  a  five-point  stencil.  However,  instead  of  using  Finite-Difference  formulae, 
the  Finite  Analytic  method  solves  exactly  the  piecewise  constant  coefficient  problems  as¬ 
sociated  with  each  cell  of  the  domain  over  which  the  governing  equation  is  defined.  The 
relations  between  neighboring  degrees  of  freedom  are,  in  effect,  matching  conditions  for 
the  solutions  continuously  defined  over  each  cell. 

The  system  of  governing  equations  described  throughout  this  section  for  use  in  model¬ 
ing  the  turbulent  reactor  could  be  solved,  after  transformation  to  streamfunction-vorticity 
form  to  eliminate  the  continuity  equation,  by  exactly  the  same  type  of  fully-implicit  mod¬ 
ified  Newton  solver  employed  in  AXIJET-L.  The  wall  function  formalism  is  slightly  dif¬ 
ferent  under  this  transformation  in  that  the  sclution  of  the  transcendental  equation  for 
the  wall  shear  stress  (or  equivalently,  for  the  friction  velocity)  involves  integral  evaluations 
when  expressed  in  terms  of  streamfunction.  It  was  originally  intended  that  AXIJET-L 
and  AXIJET-T  employ  the  same  solver,  and  simply  have  different  residual  evaluation 
routines.  However,  this  uniformity  was  not  preserved  in  the  final  product  in  favor  of  a 
more  conservative  scheme  (in  terms  of  resemblence  to  other  state  of  the  art  codes)  for 
AXIJET-T.  Apart  from  the  minor  convenience  of  handling  the  turbulence  wall  functions 
in  primitive  variables,  the  main  reasons  for  delivering  a  field-by-field  decoupled  primitive 
variable  version  were:  difficulties  with  the  k  and  c  equations  inside  of  the  Newton  loop, 
the  strongly  convection  dominated  character  of  turbulent  flows  (as  opposed  to  the  more 
elliptic  character  of  lower  Reynolds  number  laminar  flows),  which  makes  the  decoupling 
penalty  less  significant  in  this  regime,  and  the  inconvenience  of  the  special  conditionals 
required  to  handle  the  upper-bounded  fields  of  /  and  g  in  the  Jacobian  matrix  evaluation 
routine  of  the  AXIJET-L  solver.  The  nonlinear  source  terms  of  the  k  and  c  equations  are 
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not  in  and  of  themselves  worse  than  the  Arrhenius-type  source  terms  successfully  handled 
by  a  damped-modified  Newton  method  in  AXIJET-L.  However,  the  mean-flow  dissipation 
terms  (represented  by  $)  represent  a  new  spatial  coupling  between  the  mean  momentum 
and  turbulence  equations  which  is  likely  at  least  part  of  the  difficulty,  especially  since  the 
required  velocity  gradients  strongly  involve  the  “corner  points”  of  the  nine-point  discrete 
stencil.  In  practice,  Newton’s  method  was  forced  to  proceed  along  such  tiny  time-steps 
when  given  a  “reasonably”  constructable  initial  estimate  for  the  solution  iterate  that  the 
SIMPLER-type  scheme  was  superior  in  some  tests  of  non-reacting  flows.  Ultimately,  for 
large  numbers  of  species  per  gridpoint,  Newton’s  implicit  advantages  must  prevail  over 
segregated  algorithms,  and  may  yet  prevail  even  at  modest  numbers  of  species  such  as  the 
current  AXIJET-T,  but  this  is  left  to  future  research. 

One  advancement  of  the  state  of  the  art  in  primitive  variable  solvers  is  available  which 
is  not  included  in  the  current  version  of  AXIJET-T:  namely  fully  implicit  velocity-pressure 
coupling,  as  advocated  by,  t.g Vanka  [59],  We  have  successfully  tested  a  fully  implicit 
velocity-pressure  algorithm  employing  a  direct  solver  and  Newton’s  method  on  a  laminar 
non-reacting  flow  problem,  but  this  is  not  yet  incorporated  into  the  u  —  v  —  p  portion  of 
AXIJET-T.  We  do,  however,  overcome  one  of  the  prime  inefficiencies  of  the  semi-implicit 
velocity-pressure  coupling  of  the  SIMPLER  algorithm,  in  that  we  do  employ  a  direct 
(sparse)  solver  for  the  pressure  and  pressure-correction  equations. 

4.2.6.  Computational  Results 

We  show  in  this  section  results  from  three  problems  run  using  AXIJET-T  which  are 
small  enough  to  run  on  a  micro  VAX:  a  non-reacting  flow  in  the  design  configuration  of  two 
coaxial  jets,  a  reacting  flow  in  the  same  configuration,  and  a  reacting  flow  in  an  alternate 
configuration  with  one  axial  and  one  radial  inlet. 

The  former  illustrates  the  example  in  the  AXIJET-T  User’s  Guide,  a  typical  labora¬ 
tory  configuration  comprised  of  a  higher-speed  narrow  methane  fuel  jet  (0jr  =  8.23  m/s, 
rin  =  10  cm,  T  =  332  K)  and  a  low-speed  wide  oxidizer  stream  of  standard  atmospheric 
composition  ( Vq  =  0.67  m/s,  rj„  =  10  cm,  T  =  283  K).  The  inlet  conditions  were  chosen  to 
be  in  stoichiometric  balance  and  to  produce  an  exit  Reynolds  number  based  on  the  diam¬ 
eter  of  104,  which  implies  for  these  dimensions  and  inlet  conditions  a  fuel  inlet  Reynolds 
number  of  7,779  and  an  oxidizer  inlet  Reynolds  number  of  9,328.  The  exit  is  at  ZoUt  =  1.5 
m,  or  15  radii,  downstream. 

A  little  less  than  1  hr  of  MicroVAX  time  was  required  to  obtain  the  solution  depicted 
in  Figures  30  through  33  on  a  fairly  coarse  grid  (30  x  25),  using  a  more  conservative  (and 
thus  less  efficient)  than  necessary  time-step  strategy  in  which  an  initial  step  of  10-3  was 
increased  in  steps  to  10“ 1.  The  solution  is  not  yet  grid-independent  at  this  resolution,  but 
the  problem  provides  an  adequate  test  of  the  algorithmic  working  of  the  code.  The  second, 
on  the  same  grid,  is  under-resolved  like  the  first,  given  the  sharp  density  differences  of 
hydrocarbon  combustion,  but  is  included  as  a  companion  test  problem  involving  the  non- 
monotonic-with-/  state  relationship  characteristic  of  hydrocarbon-air  flames.  Its  solution 
(also  requiring  under  an  hour  of  MicroVAX  time,  under  the  same  time-stepping  strategy) 
is  represented  in  Figures  34  through  36.  The  same  exit  Reynolds  number  of  104  requires 
a  greater  inlet  mass  fluxes  due  to  the  high  viscosity  and  low  density  of  the  hot  burned 
mixture.  The  inlet  velocities  are  Vp  =  30.94  m/s  and  Vq  =  2.53  m/s,  leading  to  Reynolds 
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numbers  of  29,231  for  the  fuel  and  35,051  for  the  oxidizer. 

Reference  to  a  third  problem  is  included  briefly  in  the  report,  but  not  in  the  User’s 
Guide  because  it  represents  a  custom  application  of  AXIJET-T  to  a  very  difficult  problem 
in  the  oxidation  of  a  metallic  chloride  in  heated  air.  The  extreme  density  ratio  in  the 
problem  is  about  15.6,  which  is  roughly  twice  the  ratio  of  atmospheric  methane-air  flames, 
and  there  is  a  large  recirculation  region  tucked  out  of  the  way  of  two  cross- jets.  A  streamline 
and  temperature  plot  are  given  in  Figures  37  and  38.  The  resolution  of  the  solutions 
pictured  is  90  x  60,  which  has  been  shown  to  provide  grid-independence  for  all  functionals 
of  engineering  interest. 
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TABLE  1 


Reaction  Mechanism  Rate  Coefficients  In  The  Form  kj  =  AT^exp(- Eq/ RT). 
Units  are  moles,  cubic  centimeters,  seconds,  Kelvins  and  calories/mole. 


REACTION 

A 

P 

E 

1. 

C  H\  +  M  CH3  +  H  +  M 

1.00E+17 

0.000 

86000. 

2. 

CH\  +  02  v-  CH3  +  H02 

7.90E4-13 

0.000 

56000. 

3. 

ch4  +  h^  ch3  4-  h2 

2.20E+04 

3.000 

8750. 

4. 

CH4  4-  O  **  CH3  4-  OH 

1.60E+06 

2.360 

7400. 

5. 

CH4  +  OH  **  CH3  +  H20 

1.60E4-06 

2.100 

2460. 

6. 

CH20  4  OH  ^  HCO  +  H20 

7.53E+12 

0.000 

167. 

7. 

CH20  4-  H  **  HCO  +  H2 

3.31E4-14 

0.000 

10500. 

8. 

CH20  +  M  ^  HCO  +  H  +  M 

3.31E4-16 

0.000 

81000. 

9. 

CH20  +  O  **  HCO  4-  OH 

1.81E+13 

0.000 

3082. 

10. 

HCO  +  OH  *=*  CO  4-  H20 

5.00E+12 

0.000 

0. 

11. 

HCO  +  M*±H  +  CO  +  M 

1.60E4-14 

0.000 

14700. 

12. 

HCO  4-  H  ^  CO  +  H2 

4.00E4-13 

0.000 

0. 

13. 

HCO  4-0  **  OH  +  CO 

1.00E+13 

0.000 

0. 

14. 

HCO  +  02  ^  H02  +  CO 

3.00E4-12 

0.000 

0. 

15. 

CO  4  O  4  M  ^  C02  4"  M 

3.20E4-13 

0.000 

-4200. 

16. 

CO  A  OH  ^  C02  4-  H 

1.51E4-07 

1.300 

-758. 

17. 

C0  +  02*±  C02  4-  O 

1.60E4-13 

0.000 

41000. 

18. 

CH3  +  02  ^  CH30  4-  O 

7.00E4-12 

0.000 

25652. 

19. 

CH30  4 -M  ^  CH20  +  H  +  M 

2.40E+13 

0.000 

28812. 

20. 

CH30  +  H  ^  CH20  4-  H2 

2.00E4-13 

0.000 

0. 

21. 

CH30  4 -OH  ^  CH20  4-  H20 

1.00E+13 

0.000 

0. 

22. 

CH30  +  0*=±  CH20  4-  OH 

1.00E+13 

0.000 

0. 

23. 

CH30  4-  02  ^  CH20  +  H02 

6.30E4-10 

0.000 

2600. 

24. 

CH3  +  02  ^  CH20  +  OH 

5.20E4-13 

0.000 

34574. 

25. 

CH3  +  0^  CH20  4-  H 

6.80E4-13 

0.000 

0. 

26. 

CH3  4 -OH  **  CH20  4-  H2 

7.50E4-12 

0.000 

0. 

27. 

H02  +  CO  ^  C02  4-  OH 

5.80E4-13 

0.000 

22934. 

28. 

H2  4-  02  «=±  20  H 

1.70E+13 

0.000 

47780. 

29. 

OH  +  H2^  H20  4-  H 

1.17E4-09 

1.300 

3626. 

30. 

H  +  02^  OH  +  0 

2.20E4-14 

0.000 

16800. 

40 


TABLE  1  (continued) 

Reaction  Mechanism  Rate  Coefficients  In  The  Form  kj  =  AT^exp(—EQ/RT). 
Units  are  moles,  cubic  centimeters,  seconds,  Kelvins  and  calories /mole. 


REACTION 

A 

0 

E 

31. 

O  +  H2^  OH  +  H 

1.80E+10 

1.000 

8826. 

32. 

H  +  02  + H02  +  Ma 

2.10E+18 

-1.000 

0. 

33. 

H  +  02  +  02  *=*  H02  +  O2 

6.70E+19 

-1.420 

0. 

34. 

h  +  o2  +  n2  ho2  +  n2 

6.70E+19 

-1.420 

0. 

35. 

OH  +  H02  **  H20  +  02 

5.00E-f-13 

0.000 

1000. 

36. 

H  4-  H02  v*  2 OH 

2.50E+14 

0.000 

1900. 

37. 

0  +  ho2  ^  o2  +  oh 

4.80E+13 

0.000 

1000. 

38. 

2  OH  ^0  +  H20 

6.00E+08 

1.300 

0. 

39. 

H2  +  M^H  +  H  +  Mb 

2.23E+12 

0.500 

92600. 

40. 

o2  +  m*=*o  +  o  +  m 

1.85E+11 

0.500 

95560. 

41. 

H  +  OH  +  M  *■»  H20  +  Mc 

7.50E+23 

-2.600 

0. 

42. 

H  +  H02  ^H2  +  02 

2.50E+13 

0.000 

700. 

°  Third  body  efficiencies:  k32(H20)  =  21fc32(Ar),  £32(^2)  =  3.3A:32(Ar), 

*32(^2)  =  k32[°2)  =  0. 

b  Third  body  efficiencies:  k39(H20)  =  6k^g(Ar),  k39(H)  =  2k39[Ar),  k$9(H2)  =  3k39(Ar). 
c  Third  body  efficiency:  k4i(H20)  =  20fc4j(i4r). 
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TABLE  2 


0 . OOOOE+OO  0 . 1242E-02  0.1788E-03  1 

0 . 4721E-05  0 . 1240E-02  0.1790E-03  2 

0 . 2993E-04  0.1231E-02  0.1799E-03  3 

0 . 5602E-04  0.1222E-02  0.1809E-03  4 

0.1143E-03  0 . 1201E-02  0.1831E-03  5 

0.2398E-03  0.1162E-02  0.1875E-03  8 

0 . 4946E-03  0.1092E-02  0.1959E-03  7 

0 . 9725E-03  0.9865E-03  0.2103E-03  8 

0 . 1784E-02  0 . 8524E-03  0.2327E-03  9 

0.3027E-02  0 . 7106E-03  0.2632E-03  10 

0 . 4765E-02  0.5817E-03  0.3007E-03  11 

0 . 7030E-02  0.4754E-03  0.3432E-03  12 

0 . 8389E-02  0 . 4311E-03  0.3656E-03  13 

0.9899E-02  0.3923E-03  0.3885E-03  14 

0 . 1156E-01  0 . 3586E-03  0.4114E-03  15 

0 . 1337E-01  0 . 3295E-03  0.4341E-03  16 

0 . 1742E-01  0.2832E-03  0.4776E-03  17 

0 . 1968E-01  0 . 2650E-03  0.4977E-03  18 

0 . 2208E-01  0 . 2495E-03  0.5167E-03  19 

0 . 2464E-01  0 . 2361E-03  0.5345E-03  20 

0 . 2735E-01  0 . 2246E-03  0.5510E-03  21 

0 . 3022E-01  0 . 2147E-03  0.5662E-03  22 

0 . 3323E-01  0 . 2061E-03  0.5802E-03  23 

0 . 3640E-01  0 . 1986E-03  0.5931E-03  24 

0.3970E-01  0 . 1919E-03  0.6048E-03  25 

0.4314E-01  0 . 1861E-03  0.6156E-03  26 

0 . 4670E-01  0.1808E-03  0.6254E-03  27 

0 . 5036E-01  0 . 1761E-03  0.6344E-03  28 

0 . 5410E-01  0 . 1717E-03  0.6427E-03  29 

0 . 5866E-01  0.1672E-03  0.6514E-03  30 

0 . 6099E-01  0 . 1653E-03  0.6546E-03  31 

0 . 6337E-01  0 . 1639E-03  0.6566E-03  32 

0 . 6666E-01  0 . 1628E-03  0.6572E-03  33 

0.7014E-01  0 . 1625E-03  0.6555E-03  34 

0 . 7773E-01  0 . 1635E-03  0.6475E-03  35 

0 . 8615E-01  0 . 1654E-03  0.6371E-03  36 

0 . 9539E-01  0 . 1678E-03  0.6252E-03  37 

0 . 1055E+00  0 . 1704E-03  0.6124E-03  38 

0 . 1167E+00  0 . 1734E-03  0.5987E-03  39 

0 . 1320E+00  0 . 1775E-03  0.5805E-03  40 

0 . 1492E+00  0 . 1821E-03  0.5611E-03  41 

0 . 1684E+00  0 . 1872E-03  0.5406E-03  42 

0.1898E+00  0.1929E-03  0.5190E-03  43 

0.2137E+00  0 . 1993E-03  0.4964E-03  44 

0 . 2403E+00  0 . 2063E-03  0.4730E-03  45 

0.2699E+00  0 . 2142E-03  0.4487E-03  46 

0 . 3028E+00  0 . 2230E-03  0.4238E-03  47 

0.3393E+00  0.2329E-03  0.3984E-03  48 

0 . 3795E+00  0 . 2439E-03  0.3727E-03  49 

0 . 4236E+00  0.2563E-03  0.3470E-03  50 

0 . 4718E+00  0.2703E-03  0.3214E-03  51 

0 . 5237E+00  0 . 2862E-03  0.2962E-03  52 

0.5787E+00  0 . 3041E-03  0.2718E-03  53 

0 . 6359E+00  0 . 3246E-03  0.24S5E-03  54 

0 . 6939E+00  0 . 3480E-03  0.2265E-03  55 

0 . 7507E+00  0.3746E-03  0.2062E-03  56 

0 . 8043E+00  0.4045E-03  0 . 1877E-03  57 

0 . 8524E+00  0.4375E-03  0.1715E-03  58 

0 . 8933E+00  0 . 4724E-03  0.1575E-03  59 

0 . 9262E+00  0 . 5073E-03  0.1460E-03  60 

0 . 9511E+00  0 . 5400E-03  0.1369E-03  61 

0 . 968  9E+00  0 . 5683E-03  0.1301E-03  62 

0 . 9755E+00  0 . 5803E-03  0.1275E-03  63 

0 . 1000E+01  0 . 5886E-03  0.1246E-03  64 

-1.0  0.  0. 

F  DENSITY  VISCOSITY 


(end  of  input) 


This  file  contains  the  state  functions 
of  density  and  viscosity  to  be  read  by 
program  BETPDF.  These  data  are  from  the 
Tsuji  counter flow  quoted  from  Keyes  & 
Smooke,  J.  Comp.  Phys.,  v.73,  pp. 267-288 
(1987) .  The  card  beginning  with  "-1.0" 
is  the  sentinel  required  by  BETPDF. 
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Figure  2:  Schematic  of  an  unconfined  axisymmetric 
laminar  diffusion  flame. 
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R  IN  CM 
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Figure  3:  Experimental  (o)  and  computational  (solid 
line)  temperature  profiles  for  the  confined  methane-air 
laminar  diffusion  flame  at  a  height  of  1.2  cm  above  the 
burner  inlet. 
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MOLE  FRACTIONS 


Figure  4:  Comparison  between  measured  CJ/4  (n), 
C>2>  (°)  and  N2  (o)  profiles  and  corresponding  compu¬ 
tational  values  (solid  line)  at  a  height  of  1.2  cm  above 
the  burner  inlet. 
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MOLE  FRACTIONS 


Figure  5:  Comparison  between  measured  H2O  (n), 
CO2  (°),  CO  (o)  and  H2  (+)  profiles  and  corresponding 
computational  values  (solid  line)  at  a  height  of  1.2  cm 
above  the  burner  inlet. 
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R  IN  CM 


Figure  6:  Experimental  (o)  and  computational  (solid 
line)  temperature  profiles  for  the  confined  methane-air 
laminar  diffusion  flame  at  a  height  of  2.4  cm  above  the 
burner  inlet. 
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MOLE  FRACTIONS 


Figure  7:  Comparison  between  measured  CH4  (n), 
O2 ,  (o)  and  N2  (o)  profiles  and  corresponding  compu¬ 
tational  values  (solid  line)  at  a  height  of  2.4  cm  above 
the  burner  inlet. 
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MOLE  FRACTIONS 


Figure  8:  Comparison  between  measured  II2O  (n), 
CO2  (°).  (o)  and  IIj  (+)  profiles  and  corresponding 

computational  values  (solid  line)  at  a  height  of  2.4  cm 
above  the  burner  inlet. 
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Figure  9:  Temperature  isotherms  for  the  confined 
methane-air  laminar  diffusion  flame. 
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Mass  Fraction 


Figure  12:  Water  (Hz0)  isopleths  for  the  confined 
methane-air  laminar  diffusion  flame. 
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Figure  13:  Carbon  Monoxide  (CO)  isoplcths  for  the 
confined  methane-air  laminar  diffusion  flame. 
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Figure  14:  Molecular  hydrogen  (H2)  isopleths  for  the 
confined  methane-air  laminar  diffusion  flame. 
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Figure  15:  Carbon  Dioxide  (C02)  isopleths  for  the 
confined  methane-air  laminar  diffusion  flame. 
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Figure  17:  Hydroxyl  radical  (OH)  isopleths  for  \ 
confined  methane-air  laminar  diffusion  flame 
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Figure  18:  Formaldehyde  (CH20)  isopleths  for  the 
confined  mcthanc-air  laminar  diffusion  flame. 
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Figure  19:  Stream  function  isopleths  for  the  confined 
methane-air  laminar  diffusion  flame. 
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Figure  20:  Temperature  isotherms  for  the  unconfined 
methane-air  laminar  diffusion  flame. 


62 


Streamlines 


0.091 


0.09 


0.00 


0.07 


0.06 


0.05 


0.04 


0.03 


0.02 


0.01 


0.00 


20.0 


Figure  21:  Stream  function  isopleths  for  the  unconfined 
methane-air  laminar  diffusion  flame. 
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Figure  22:  Vorticity  isopleths  for  the  unconfined 
mcthanc-air  laminar  diffusion  flame. 
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Figure  25:  Fuel  mass  fraction  (peaking  at  /  =  1),  oxy¬ 
gen  mass  fraction  (peaking  at  /  =  0),  and  temperature 
(peaking  near  /  =  fst)  for  the  Tsuji  counter  flow. 


FUEL 


Figure  26:  Schematic  of  confined  coaxial  configuration 
solved  by  AXIJET-T. 
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VISCOUS  SUBLAYER 


Figure  27:  versus  u'*’  on  a  linear-log  plot,  illustrating 

the  log  law. 
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Figure  28:  fc+  versus  y+,  illustrating  the  wall  law  for 

k. 
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Grid  Distribution 
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Figure  JO:  Grid  distribution  for  the  cold  mixing  prob¬ 
lem. 
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Grid:  30x25;  Unit:  inch 
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H.-0.900E+01  15 : 0  .  lOOE-t  02 
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Figure  SI:  Streamfunction 
mixing  problem 


contours  for  the  methane-air 
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Grid: 


( inch) 


Conserved  Scalar  Contours 
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Figure  32:  Mixture  fraction  contours  for  the  methane- 
air  mixing  problem. 
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Grid:  30x25; 


( inch) 


Turbulent  Dissipation  Rate  Contours 
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Figure  33:  Turbulent  kinetic  energy  dissipation  rate 
contours  for  the  methane-air  mixing  problem. 
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Figure  34:  Streamfunction  contours  for  the  methane-air 
flame 
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Grid:  30x25: 


( inch) 
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Figure  35:  Mixture  fraction  contours  for  the  methane- 
air  flame. 
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Grid:  30x22; 


( inch) 


Turbulent  Dissipation  Rate  Contours 
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Figure  30:  Turbulent  kinetic  energy  dissipation  rate 
contours  for  the  methane-air  flame. 
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Figure  37:  Streamline  contours  for  a  custom  AXIJET-T 
application  with  cross  jets,  and  large  density  ratio,  on 
a  90  by  60  grid. 
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Grid:  90x00; 


( inch) 
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Figure  38:  Temperature  contoure  for  a  custom  AXI.IET- 
T  application  with  cross  Jets,  and  large  density  ratio,  on 
a  90  by  00  grid. 
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Figure  SO:  Sample  beta  function  probability  distribu¬ 
tions  at  small  mixture  fraction  covariance. 
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82 


5.1.  Introduction 

In  practical  calculations  of  turbulent  reacting  flows,  Favre-averaged  values  of  the  mix¬ 
ture  properties  are  often  evaluated  through  the  introduction  of  a  probability  density  func¬ 
tion  (pdf)  for  a  single  conserved  scalar,  /.  The  Favre  average  is  a  density-weighted  time 
average  of  a  fluctuating  quantity,  and  is  customarily  denoted  by  means  of  a  tilde  to  dis¬ 
tinguish  it  from  the  overbar  employed  in  ordinary  time-averaging.  For  a  fluctuating  field 

4>, 


i  =  p4>!p- 


(5.1) 


The  pdf  P(/,x)  when  multiplied  by  a  differential  d /  gives  the  probability  finding  the 
instantaneous  value  of  the  conserved  scalar  in  a  range  d /  about  /  at  the  point  x.  The 
total  range  of  the  conserved  scalar  is  0  <  /  <  1.  Thus,  the  Favre  average  of  any  quantity 
<f>  depending  solely  on  /,  such  as  a  species  concentration  under  an  equilibrium  kinetics 
assumption,  for  instance,  is  given  at  a  point  x  by 


f\(f)P(f,x)df.  (5.2) 

JO 

In  this  expression  is  the  equilibrium  value  of  the  quantity  <f>  corresponding  to  the 
instantaneous  /. 

The  pdf  P(f,x)  is  assumed  to  embody  all  aspects  of  the  local  turbulent-chemistry 
interaction.  Under  the  chemical  equilibrium  hypothesis,  as  well  as  under  other  hypotheses 
not  illustrated  herein,  a  useful  conserved  scalar  is  the  mixture  fraction,  defined  as 


Z  —  Zf 

Zq  -  Zp' 


(5.3) 


where  Z  is  any  conserved  function  of  the  reaction  (such  as  a  stoichiometrically  weighted 
linear  combination  of  fuel  and  oxidant),  and  where  subscripts  F  and  O  denote  its  values 
in  the  fuel  and  oxidant  streams,  respectively.  The  form  of  an  incomplete  beta  function  is 
often  adopted  for  the  pdf.  Its  form  was  given  in  equations  (4.92)  and  (4.93). 

The  Favre  averages  /  and  g  at  the  point  x  are  obtained  by  solving  modeled  partial 
differential  equations  of  transport  type  as  part  of  the  main  computational  procedure,  which 
often  includes  a  k  —  e  turbulence  model.  The  resulting  formalism  is  sometimes  referred 
to  as  the  k  -  c  -  g  model.  A  thermochemical  package  is  needed  to  provide  <f>(f)  for  each 
species  concentration,  viscosity,  temperature,  etc.  Substitution  of  (4.92)  into  (5.2)  gives 


83 


(5.4) 


i = wj)  L' -  f)>~' 4f' 

where  the  complete  beta  function  is  defined  by  the  definite  integral 

0(«,6)=  (  /"'‘(I -/)*■'  H- 

Jo 

The  formula  for  the  gas  density  is  of  a  slightly  different  type  since  the  time-averaged  density 
is  included  in  the  Favre-averaged  definition  of  all  other  fields.  The  appropriate  relationship 
using  the  same  pdf  (4.92)  is 

w~1  =  j<h ifam-'r-'H-it-'i/.  <»■*> 

Equations  (5.4)  and  (5.5)  are  both  of  the  type 

|  Jd^t*)**-1^-*)*-1** 

/0V-i( 

Since  the  variance  g  is  constrained  to  He  in  the  interval  0<g  <  f  (l  —  f),  a  and  6  are 
both  greater  than  zero  from  (4.93),  and  thus  the  numerator  and  denominator  of  (5.6) 
are  integrable  for  all  physically  reasonable  The  results  of  the  integration  are  best 
tabulated  once  for  repeated  “lookup”  during  the  governing  equations  solution  phase.  A 
few  considerations  appropriate  to  this  tabulation  are: 

•  When  a  and  b  attain  large  values,  numerical  quadrature  of  the  integrals  in  (5.6)  can 
lead  to  underflow,  since  ^  is  bounded,  and  at  least  one  factor  of  each  integrand  tends 
to  zero  in  this  limit.  Even  in  the  absence  of  underflow,  a  modestly  large  absolute  error 
tolerance  in  the  quadratures  could  result  in  an  unacceptably  inaccurate  <f>  being  formed 
from  the  ratio. 

•  When  a  or  6  is  less  than  one,  the  integrals  become  improper,  since  their  integrands 
tend  to  infinity  at  at  least  one  endpoint. 

•  An  efficient  arrangement  of  the  table  is  needed  to  minimize  the  CPU  time  spent  search¬ 
ing  the  table  during  the  main  computation. 

We  recommend  in  the  next  section  an  efficient  and  robust  procedure  for  dealing  with 
these  potential  problems. 
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5.2.  Numerical  Evaluation  of  Turbulent  State  Relations 

In  most  cases  of  practical  interest,  the  relationship  <£(x)  between  the  mixture  property 
and  the  instantaneous  conserved  scalar  is  given  in  the  form  of  tabulated  data,  that  is,  4>{x) 
is  known  only  at  a  set  of  n  +  1  points  {x,},  0  =  xq  <  x^  <  •••  <  xn-\  <  xn  =  1,  which 
partition  the  unit  interval  into  n  subintervals.  On  each  subinterval,  4>{x)  may  conveniently 
be  approximated  by  a  polynomial  of  degree  k  whose  coefficients  depend  on  the  values  ^(xt) 
at  some  neighboring  set  of  k+1  indices.  For  piecewise  linear  interpolation  in  Xj_j  <  x  <  xt-, 


4>{x)  =  <t>{xi_x)  + 


g  -  g»- 1 

xt  ~  xi— 1 


(<£(xt)  -  ^(X,-!)), 


or 


<£(x)  =  e,x  +  <^,  (5.7) 

Only  the  piecewise  linear  case  will  be  explicitly  treated  in  the  sequel;  however,  generaliza¬ 
tion  to  any  finite-order  polynomial  can  be  straightforwardly  accommodated  in  the  following 
development  by  the  inclusion  of  additional  terms  whose  integrands  contain  higher  powers 
of  x.  No  cases  requiring  special  treatment  occur  for  powers  beyond  the  zeroth  and  first. 

Employing  the  interpolation  (5.7)  in  the  numerator  of  (5.6)  gives 


f  <t>(x)xa  Hi-*)6  ldx  =  ^2  Ci  f 
Jo  ,=1  Jxi- 1 


xa(l  —  x)b  1  dx  +  d,- 


f 

JXi-l 


x)b~l  dx 


=  j^k(0(«  +  1  ,&.*»•)  -/?(o  +  l,6,x,_!))  +  d,(/?(a,6,x,)  -/?(«,  6,  x^))] , 


t  =  l 


where 


/?(a, 6, x)  =  f  xa  *(1  -  x)b  1  dx 
JO 

is  the  unnormalized  incomplete  beta  function.  (Note  that  in  this  triple-argument  notation 
0[a,b,  1)  is  simply  the  complete  beta  function.)  The  denominator  of  (5.6)  can  be  handled 
in  exactly  the  same  manner,  with  Cj  =  0  and  d,  =  1  for  i  =  1,  •  •  • ,  n. 


Equation  (5.6)  can  also  be  written  as 


(i  -  *)*-> 


dl  =  53  Mi  +  «UI  • 
*=1 


(5.8) 
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where 


lx  =  /(a  +  1,6, xt-)  -  /(a  +  1,  *,**_!),  Jt  =  /(a,6,x<)  -  /(a.M^), 
where,  in  turn, 

/(a’M  =  WJ)  /„  X°~'{1  ~ 1)1-1  iz  ,59) 

is  the  standard  notation  for  the  normalized  incomplete  beta  function. 

Since  the  It  and  are  independent  of  ^(x),  they  may  be  calculated  once  at  each 
(a,  6)  required  in  the  state  function  tables.  The  state  function  tables  for  each  field  are  then 
produced  by  summing  over  the  appropriate  set  of  e,  and  d 

The  quadratures  can  be  done  by  either  of  two  convenient  methods.  If  accurate  sub¬ 
routines  for  the  evaluation  of  (5.9)  are  available,  they  may  be  employed  directly.  If  not, 
adaptive  quadrature  using  a  Richardson-extrapolated  Simpson’s  formula  on  successively 
refined  meshes  will  work  in  subintervals  2  through  n  —  1.  In  the  BETPDF  program,  this 
is  embodied  in  the  FORTRAN  routine  QSIMP. 

As  noted  in  the  introduction  section,  there  are  two  regimes  of  the  exponential  param¬ 
eters  of  the  beta  function  wnich  demand  special  numerical  treatment. 

•  The  case  of  a  and  6  sufficiently  large  to  cause  underflow  occurs  only  if  g  is  very  small. 

Physically,  this  implies  that  the  variances  of  the  fluctuating  state  variables  about  their 
means  is  negligible.  In  this  limit,  the  pdf  P{f)  may  be  represented  by  a  Dirac  delta 
function  centered  at  /.  Hence,  ^  directly,  and  numerical  integration  is  unnec¬ 

essary. 

•  If  either  a  or  b  is  less  than  unity,  the  left-  or  right-most  subintervals,  respectively, 
require  a  preliminary  integration  by  parts  to  remove  the  negative  power  of  xa~l  or 
(1  —  x)6-1,  respectively.  Thus,  one  obtains  the  special  formulae: 

J\  =  -(xi)a(l  -  xi)i_1  +  - — ^  /  xfl(l  -  x)b~2  dx,  when  0  <  a  <  1, 

a  a  J  o  k 

Irx  =  x(xn— i)“(l  -  xn— i)6  +  ~  f  xa-1(l  -  x)b  dx,  when0<6<l,  * 

°  b  J*n-\ 

J n  =  T(xn_i)0-1(l  -  xn_i)ft  +  — r—  /  x0-2(l  -  x)ft  dx  when  0  <  b  <  1. 

°  6  J*n- 1 

(Note  that  no  special  formula  is  needed  for  Jj  regardless  of  the  magnitude  of  a.) 
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•  As  further  noted  in  the  introduction,  the  state  function  routines  will  be  called  repet¬ 
itively  at  each  grid  point  at  each  iteration  during  the  main  computational  procedure 
which  solves  the  governing  Favre-averaged  partial  differential  equations.  Thus,  the 
table  searches  need  to  be  very  efficient  in  terms  of  operation  count.  The  domain 
0<  /  <  1,  0  <  $  <  /(I  —  /)  has  a  parabolic  bounding  arc.  This  implies  that  a  tabula¬ 
tion  of  state  data  over  a  uniform  tensor  product  mesh  in  the  circumscribing  rectangle  of 
(f,g)- space,  implemented  as  a  two-dimensional  FORTRAN  array,  is  wasteful  of  points. 
Any  other  type  of  mesh  requires  an  iterative  search  in  g  followed  by  a  bilinear  inter¬ 
polation  procedure  in  a  nonrectangular  region.  However,  the  transformation  to  (/,  q) 
coordinates,  where 

" =  /(i  -  f) 

and  0</<l,0<q<l  makes  the  state  space  a  square.  In  terms  of  q , 

a  =  /( i  -  1),  b  =  (1  -  /)(i  -  1). 

n  n 

Tabulation  of  the  state  variables  at  (/,,  qy),  where  fa  =  i/Nf,  i  =  0,  •  •  • ,  Nf  and  where 
qy  =  j/Nn,  j  =  0, provides  that  ^(/,q)  may  be  bilinearly  interpolated  in 
exactly  eight  scalar  multiply-adds  and  one  scalar  division. 

S.3.  Example  and  Discussion 

The  calculation  of  the  Favre-averaged  viscosity  and  density  in  terms  of  /  and  q  is 
presented  as  an  example.  The  equilibrium  data  for  density  and  viscosity  is  listed  in  Table  II. 
Figures  23  and  24  display  three-dimensional  views  of  pdf-averaged  density  and  viscosity 
on  the  (/,q)  plane,  respectively.  (/  =  0  corresponds  to  the  oxidizer  stream  and  /  =  1  to 
the  fuel.) 

The  striking  features  of  these  figures  are  the  non-monotic  relationships  between  the 
state  variables  and  the  mixture  fraction  at  small  q  and  the  asymptotic  approach  of  the 
viscosity  curve  to  a  straight  line  between  the  fuel  and  oxidant  reservoir  values  when  q  — »  1. 
The  first  feature  is  not  necessarily  characteristic  of  all  fuel-oxidizer  combinations  and  inlet 
temperatures.  (Hydrocarbon-air  mixtures  are  notorious  for  the  nearly  vertical  tangents 
of  the  density  at  small  values  of  /,  where  /  passes  through  its  stoichiometric  value,  and 
this  fact  needs  to  be  taken  into  account  in  adaptively  refining  the  resolution  of  discrete 
methods  in  such  regions.  Failure  to  do  so  will  nullify  the  accuracy  of  any  calculation  in 
which  density  differences  play  a  role.)  The  second  (highly  appropriate)  feature  is  generic 
to  all  beta-function  generated  models,  can  be  verified  mathematically  as  follows. 

In  the  limit  of  q  — ►  1,  we  may  write  q  =  1  —  c  and  consider  instead  the  limit  e  — ►  0 
(through  positive  values).  From  a  Taylor  expansion  with  second-  and  higher-order  terms 
neglected,  we  then  have  osj  q/  and  6  «  q(l  -  /).  From  the  alternative  integral  definition 
of  the  complete  beta  function, 


0(a,6)  =  f\  1  +  V1  +  *6_1)  d*. 

JO 

we  easily  obtain  to  leading  order  that  0(atb)  «  1/e  in  this  limit.  Thus,  the  denominator 
of  (5.6)  becomes  unbounded.  Let  the  numerator  be  written 


* 


rS  rl—6  f 1 

/  ^(x)xa-1(l  -  x)b~l  dx  +  /  ^(x)x0-1(l  -  x)6_1  dx  +  /  <£(x)xa-1(l  -  x)fc_1  dx, 

Jo  Js  Jl-6 

where  6  is  some  fixed  positive  constant  less  them  unity  which  is  chosen  sufficiently  small 
that  can  be  represented  to  adequate  accuracy  by  simply  the  constant  term  in  its 
Taylor  series  about  the  extreme  endpoints  in  the  first  and  third  integrals.  (Since  the  rest 
of  each  integrand  has  delta  function-like  behavior  at  each  endpoint  as  e  — ►  0,  this  is  always 
possible.)  Then  the  first  and  third  terms  of  the  numerator  become,  respectively, 


6(f 

4(0)-? +  *W 

e/ 


6<(W) 
<(1  -  /) 


The  middle  term  of  the  numerator  can  be  bounded  independently  of  e  by 


max 


1  -  26 
6 2 


Multiplying  both  numerator  and  denominator  by  e  and  letting  e  -♦  0  gives  the  result  that 


4>  -*  MO)  +  (1  -  /)<£(!), 


as  was  to  be  shown. 

Evaluation  of  Favre-averaged  properties  of  a  turbulent  reacting  flowfield  using  a  beta 
function  formalism  is  mathematically  straightforward  but  requires  some  care  in  implemen-  » 

tation  in  finite  precision  arithmetic,  especially  considering  the  strong  nonlinear  couplings 
in  which  the  state  variables  are  involved  in  the  system  of  governing  equations.  In  the 
technique  described  above  for  the  accurate  evaluation  of  these  properties,  the  improper 
integrals  occuring  in  the  beta  function  integrals  are  removed  by  partial  integration  before 
numerical  quadrature  is  employed.  All  of  the  mixture  properties  can  be  evaluated  from 
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one  set  of  evaluations  of  the  beta-function  integrals.  The  searching  of  the  data  table  is  sim¬ 
plified  by  means  of  a  transformation  which  makes  the  domain  of  the  independent  variables 
a  square. 


j 
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